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Preface 


The  purpose  of  this  research  was  to  deveioo  and  validate  a  fallout 
prediction  method  using  variable  winds  for  particle  transport  calcula¬ 
tions.  The  new  method  uses  National  Meteorological  Center  (NMC)  spectral 
coefficients  to  compute  wind  vectors  along  the  space  and  time  varying 
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Abstract 

A  new  method  has  been  developed  to  incorporate  variable  winds  into 
fallout  transport  calculations.  The  method  uses  spectral  coefficients 
derived  by  the  National  Meteorological  Center.  Wind  vector  components  are 
computed  with  the  coefficients  along  the  trajectories  of  falling  particles. 
Spectral  winds  are  used  in  the  two-step  method  to  compute  dose  rate  on  the 
ground,  downwind  of  a  nuclear  cloud.  First,  the  hotline  is  located  by  com¬ 
puting  trajectories  of  particles  from  an  initial,  stabilized  cloud,  through 
spectral  winds,  to  the  ground.  The  connection  of  particle  landing  points 
is  the  hotline.  Second,  dose  rate  on  and  around  the  hotline  is  computed  by 
analytically  smearing  the  falling  cloud's  activity  along  the  ground.  The 
feasibility  of  using  spectral  winds  for  fallout  particle  transport  was  vali 
dated  by  computing  Mount  St.  Helens  ashfall  locations  and  comparing  calcula 
tions  to  fallout  data.  In  addition,  an  ashfall  equation  was  derived  for 
computing  volcanic  ash  mass/area  on  the  ground.  Ashfall  data  and  the  ash¬ 
fall  equation  were  used  to  back-calculate  an  aggregated  particle  size  dis¬ 
tribution  for  the  Mount  St.  Helens  eruption  cloud.  Further  validation  was 
performed  by  comparing  computed  and  actual  trajectories  of  a  high  explosive 
dust  cloud  (DIRECT  COURSE) .  Using  an  error  propagation  formula,  it  was 
determined  that  uncertainties  in  spectral  wind  components  produce  less  than 
four  percent  of  the  total  dose  rate  varianceT^In  summary,  this  research 
demonstrated  the  feasibility  of  using  spectral  coefficients  for  fallout 
transport  calculations,  developed  a  two-step  smearing  model  to  treat  vari¬ 
able  winds,  and  showed  that  uncertainties  in  spectral  winds  do  not  contri¬ 
bute  significantly  to  the  error  in  computed  dose  rate. 


DEVELOPMENT  AND  VALIDATION  OF  A  NEW  FALLOUT 


TRANSPORT  METHOD  USING  VARIABLE  SPECTRAL  WINDS 


I .  Introduction 


Background 

An  atmospheric  nuclear  burst  creates  a  cloud  of  vaporized  radio¬ 
active  particles  composed  of  fission  fragments,  decay  products,  unfis¬ 
sioned  bomb  fuel  and  neutron- induced  debris.  The  radioactive  particles 
rise  with  the  buoyant  nuclear  cloud,  condensing  on  microscopic  nuclei 
or  on  the  surfaces  of  other  particles.  As  condensed  particles  grow  and 
the  cloud  cools,  gravitational  forces  eventually  exceed  the  forces  in¬ 
duced  by  updraft,  and  the  radioactive  particles  fall  from  the  cloud. 
During  their  fall,  particles  are  transported  by  ambient  winds.  As  a 
result,  radioactive  fallout  can  land  far  from  the  point  of  weapon  deto¬ 
nation.  Fallout's  widespread  distribution  coupled  with  its  long-term 
nature  make  fallout  transport  modeling  an  important  part  of  nuclear 
weapon  effects  studies. 

Fallout  transport  studies  have  both  military  and  civilian  applica¬ 
tions.  Militarily,  fallout  may  be  an  important  consideration  in  studies 
of  force  recovery  and  organization  after  a  strategic  nuclear  attack  or 
during  a  protracted,  limited  nuclear  war  (69) (82) .  Fallout  models  are 
also  used  to  estimate  population  exposure  from  past  atmospheric  nuclear 
weapons  tests  (8).  In  civilian  applications,  fallout  models  are  applied 
to  study  population  survival  and  economic  recovery  prospects  following  a 
nuclear  attack  (30) (42) (56) (60) (66) (76) .  Currently,  fallout  transport 
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models  are  even  used  to  estimate  potential  climatological  effects  of 
global  nuclear  fallout  and  dust,  the  "nuclear  winter"  (105). 

Fallout  Models 

Fallout  studies  use  both  numerical  and  analytic  methods  to  predict 
particle  motion  in  the  atmosphere  and  dose  rates  on  the  ground.  The 
standard  numerical  fallout  code  for  military  applications  is  DELFIC,  the 
Defense  Land  Fallout  Interpretative  Code  (79).  DELFIC  computes  cloud 
rise,  growth,  stabilization  and  transport,  with  detailed  simulation  of 
radioactive  particle  formation.  Particle  fall  is  computed  in  a  wind 
field  that  is  (optionally)  spatially  constant,  spatially  constant  with 
temporal  updates  or  spatially  interpolated  within  a  finite  grid  of  wind 
profiles.  DELFIC  is  frequently  used  as  a  computational  benchmark  for 
simpler  analytic  models,  such  as  the  AFIT  (Air  Force  Institute  of  Tech¬ 
nology)  and  WSEG  (Weapon  Systems  Evaluation  Group)  codes  (17) (19) (78) 

(85) (86) . 

AFIT  and  WSEG  simulate  the  fallout  process  by  smearing  (depositing) 
the  falling  cloud's  radioactivity  on  the  ground  unidirectionally  downwind. 
The  analytic  codes  assume  a  constant  wind  vector  to  transport  the  cloud 
(49) (65) .  The  constant  wind  assumption  is  a  significant  source  of  error, 
particularly  when  winds  are  not  nearly  constant  in  time  and  space. 

Particle  Transport  in  Spectral  Winds 

This  report  presents  a  new  way  to  compute  particle  transport  using 
variable  winds.  The  new  method  uses  spectral  coefficients  currently 
generated  by  the  National  Meteorological  Center  (NMC)  and  soon  to  be 
generated  by  the  Air  Force  Global  Weather  Central  (AFGWC) (103) .  Spectral 
coefficients  are  derived  from  global  wind  data  collected  twice  daily. 
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nearly  simultaneously.  By  fitting  the  wind  data  with  a  truncated  series 
of  spherical  harmonics,  NMC  reduces  the  data  to  a  set  of  complex  spectral 
coefficients.  The  coefficients  are  used  in  the  spherical  harmonics  expan¬ 
sions  to  compute  wind  speeds  (spectral  winds)  anywhere  in  the  atmosphere. 

With  spectral  winds,  fallout  transport  is  modeled  by  computing  wind 
speeds  from  coefficients  at  discrete  positions  and  times  along  the  trajec¬ 
tories  of  falling  particles.  Horizontal  spatial  variability  is  naturally 
simulated  by  the  harmonic  polynomials.  There  are  coefficient  sets  for 
each  of  twelve  different  altitudes,  so  winds  (or  coefficients)  must  be 
vertically  interpolated  to  obtain  wind  speeds  between  spectral  heights. 
Each  coefficient  set  fits  the  polynomials  to  global  wind  data  at  a  fixed 
time,  so  winds  (or  coefficients)  can  be  temporally  interpolated  to  esti¬ 
mate  wind  speeds  at  times  when  data  was  not  collected. 

Using  spectral  winds  to  compute  fallout  particle  trajectories  will 
vastly  improve  the  accuracy  of  fallout  prediction  codes.  Spectral  winds 
represent  real  atmospheric  wind  data  over  the  entire  globe.  Particle 
trajectories  can  be  computed  with  continuous,  four  dimensional  varia¬ 
bility. 

Analytic  Fallout  Code  with  Variable  Winds 

Spectral  winds  can  be  used  in  any  numerical  or  analytic  fallout 
code.  The  spectral  wind  equations  are  ideally  suited  for  calculations 
of  wind-driven  particle  motion  in  a  Lagrangian  reference  frame.  A  numer¬ 
ical  fallout  code  with  spectral  wind  transport  would  be  an  ideal  fallout 
prediction  sys  .em,  coupling  a  state-of-the-art  physics  model  with  global, 
variable  winds.  However,  DELFIC,  the  standard  numerical  fallout  code  is 
so  large  and  slow  that  it  is  impractical  to  use  for  most  fallout  studies. 


especially  when  many  hypothetical  bursts  are  being  considered.  Analytic 
codes,  such  as  AFIT  and  WSEG,  are  used  for  most  fallout  studies  because 
they  are  economical,  fast-running  and  benchmarked  against  DELFIC  and 
actual  fallout  data  (19) (69) (81) (84) .  In  this  study,  spectral  wind 
transport  is  incorporated  into  an  analytic  fallout  code. 

The  AFIT  and  ^SEG  analytic  fallout  codes  were  designed  for  simpli¬ 
fied  cloud  transport  with  a  single  wind  vector.  A  new  analytic  method 
was  developed  to  compute  fallout  dose  rates  with  variable  winds  (53). 

The  new  method  is  a  two  step  approach.  First,  the  radioactive  hotline 
is  located  by  tracking  the  motion  of  particles  from  an  initial,  stabi¬ 
lized  cloud,  through  spectral  winds,  to  the  ground.  The  hotline  is  the 
locus  of  peak  activity  downwind  of  the  burst.  Second,  dose  rate  is  cal¬ 
culated  with  a  solution  to  the  two-directional  dose  rate  integral 
equation: 


D  =  k  Y  ff  /  f(x,y,t)  g(t)  dt  (1.1) 

o 


where 

D  =  dose  rate  at  (x,y)  at  one  hour  after  burst 

note:  D  is  called  dose  rate  in  this  report  to  be  consistent 
with  most  fallout  references.  Actually,  5  is  the  dose  rate 
to  air,  or  exposure  rate,  as  defined  by  the  International 
Commission  on  Radiation  Units. 

k  =  source  normalization  constant,  typically  2350  (Roentgens /hour) 
(square  miles/kiloton) 

Y  *  weapon  yield  (kilotons) 

ff  =  fission  fraction 

f(x,y,t)  3  spatial  distribution  function 

g(t)  =  activity  arrival  rate  function 


With  the  two  step  method,  cloud  radioactivity  is  analytically  smeared 
along  the  ground  on  and  around  a  curved  hotline  defined  by  explicit  par¬ 
ticle  fall  calculations  through  variable,  spectral  winds. 
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Organization  of  Report 


Chapter  II  explains  how  the  NMC  derives  spectral  coefficients,  and 
how  wind  vector  components  are  computed  from  the  coefficients.  Chapter 
III  describes  the  hotline  locator  model  that  uses  spectral  winds  to  trans¬ 
port  falling  particles  from  the  initial  cloud  to  the  ground.  A  variable 
wind  fallout  code  is  presented  in  Chapter  IV,  with  an  algorithm  to  compute 
dose  rate  at  arbitrary  field  points.  Chapter  V  shows  results  of  two  vali¬ 
dation  studies,  using  data  from  a  Mount  St.  Helens  ash  cloud  and  from  a 
cloud  produced  by  a  large  high  explosive  (not  nuclear)  atmospheric  burst. 
An  error  analysis  is  presented  in  Chapter  VI,  quantifying  the  effects  of 
uncertainties  in  wind  speed  (and  other  independent  variables)  on  dose 
rate  at  a  point. 

Appendix  A  shows  a  sample  of  spectral  coefficients  for  the  1000 
millibar  height  at  00UT  (Universal  Time)  16  Jan  82.  Appendix  B  contains 
the  derivation  of  the  equation  for  computing  spectral  winds.  Appendix  C 
illustrates  how  different  nuclear  models  position  the  stabilized  cloud 
in  the  atmosphere.  Polynomials  for  computing  the  Laurent  series  coeffi¬ 
cients  (used  to  determine  particle  sizes  arriving  downwind)  are  in  Appen¬ 
dix  D.  Appendix  E  contains  an  integration  of  cloud  activity  in  the  cross- 
wind  direction.  The  algorithm  for  off-axis  dose  rate  determination  is 
in  Appendix  F.  Appendix  G  explains  how  particle  size  distributions  are 
fitted  with  trimodal  log-normal  distributions,  using  a  constrained  mini¬ 
mization  search  in  the  nine  parameter  hyperspace  defined  by  the  sum  of 
three  log-normal  distribution  functions. 
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II.  Spectral  Wind  Calculations 


In  this  report,  spectral  winds  are  the  wind  vector  components 
obtained  from  NMC  spectral  coefficients.  This  chapter  describes  the 
source  of  the  coefficients  and  the  equations  for  computing  wind  vectors. 

NMC  Spectral  Coefficients 

Numerical  models  of  atmospheric  motion  were  developed  by  the  NMC 
and  others  to  study  the  general  circulation  of  the  atmosphere  and  to 
generate  forecasts.  Atmojpheric  motion  is  modeled  by  NMC  with  a  set 
of  non-linear  partial  differential  equations  that  describe  the  temporal 
variation  of  space-dependent  atmospheric  variables  (37) (71) .  The  par¬ 
tial  differential  equations  are  converted  to  ordinary  differential  equa¬ 
tions,  then  solved  by  finite  difference  techniques  (15) (32) (34) (83) . 
Space-dependent  atmospheric  variables  in  the  partial  differential  equa¬ 
tions  are  replaced  by  spherical  harmonic  expansions  of  the  form  (89) (93) 
(95): 

J  | al+j 

D  =  I  l  dJ  Pj(sin*)  EXP(iU)  (2.1) 

&=-J  n=|2,| 

where 

D  =  space  (and  time;  dependent  atmospheric  variable 
0 

D*  =  complex  expansion  (spectral)  coefficient 

P^  =  associated  Legendre  polynomial 

<j>  =  latitude 

X  =  longitude 

l  =  zonal  (latitudinal)  index 
n  =  ordinal  (longitudinal)  index 
J  =  truncation  limit 


Substituting  a  series  like  Eq  (2.1)  for  each  space- dependent  variable 
transforms  the  partial  differential  equations  into  a  set  of  ordinary 
differential  equations,  with  time  as  the  independent  variable.  The 
truncation  limit,  J  ,  is  selected  by  NMC  to  be  high  enough  for  accurate 
forecasts  that  are  independent  of  the  truncation  limit  (87). 

As  defined  in  the  summations,  allowable  values  of  %  and  n  are 
bounded  within  a  rhomboid  on  a  plot  of  £  vs  n.  Rhomboidal  truncation 
is  used  primarily  because  it  gives  the  same  number  of  meridional  com¬ 
ponents  for  each  zonal  index  (wave  number)  (35) (94) . 

Time  integration  of  the  differential  equations  requires  initial 
values  of  the  expansion  coefficients,  D^.  The  initial  coefficients 
are  computed  by  NMC  using  observed  global  atmospheric  data  and  the 
normalization  integral  for  spherical  harmonics  (95) .  The  initial 
spectral  coefficients  fit  the  truncated  spherical  harmonics  to  observed 
global  wind  data  that  is  collected  twice  daily.  Coefficients  can  be 
used  to  compute  the  observed  wind  data  and  to  interpolate  the  data  on 
a  spherical  earth. 

A  set  of  spectral  coefficients  is  derived  lor  each  wind  component, 
u  and  v,  at  each  of  twelve  heights  in  the  atmosphere  at  a  fixed  time. 
Figure  II-l  shows  the  spectral  heights.  Individual  sets  of  coefficients 
are  not  empirically  or  analytically  related  in  the  NMC  spectral  models 
of  atmospheric  motion.  However,  the  data  used  to  derive  the  coefficients 
comes  from  continuous  vertical  soundings,  so  spectral  winds  will  have  an 
inherent  inter-altitude  correlation  at  a  fixed  time.  The  NMC  has  demon¬ 
strated  that  the  spectral  wind  structure  is  vertically  consistent  (70) . 
Therefore,  particle  trajectories,  passing  through  spectral  winds  at 
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various  altitudes,  should  be  smooth  and  continuous. 

Appendix  A  shows  some  spectral  coefficients  for  the  1000  millibar 
height  at  00UT  (Universal  Time)  on  16  Jan  82.  They  are  for  the  u-com- 
ponent  of  wind;  63  complex  pairs  are  shown.  A  complete  set  of  coeffi¬ 
cients  for  all  twelve  atmospheric  heights,  at  one  time,  consists  of 
23808  complex  pairs.  Spectral  coefficients  are  routinely  archived  by 
NMC  at  the  National  Climatic  Center,  Asheville,  North  Carolina. 


Wind  Calculation 


With  the  spectral  coefficients  for  a  specific  date,  time  and  pres¬ 
sure  level,  wind  vector  components  can  be  computed  anywhere  in  the  atmo¬ 
sphere  using  the  spherical  harmonics  expansion  (95) (112) : 


J  |a|+J+1 

(U,V)  =  (u,v)cos(j>  =  l  l  (Uj,Vj)  Pj(sin+)  EXP(iAX)  (2.2) 

i-~J  n=|A| 

where 

(U,V)  =  pseudo-wind  components,  either  U  or  V 
(u,v)  =  spectral  wind  components,  either  u  or  v 
u  =  west  to  east 
v  =  south  to  north 

(U^,V*j)  =  spectral  coefficients,  either  U^  or  V* 


By  rearranging  the  summation  and  using  identities  for  spherical 
harmonics,  the  spectral  equations  become: 

(U,V)  =  (u,v)  cos<f> 


J  A+J+l 

=  Re  l  A  l  (Uj,vj)  P^(sin<j))  (cos AX  +  i  sinAX)  (2.3) 

A-0  n=A 
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where 

A  =  1  if  1  =  0 

A  *  2  otherwise 

See  Appendix  B  for  the  stepwise  derivation  of  Eq  (2.3)  from  Eq  (2.2). 
Appendix  B  also  shows  the  recursion  relation  for  the  associated  Legendre 
polynomials,  and  the  computational  algorithm  for  spectral  winds. 

Figures  II-2  through  11-13  show  fields  of  wind  vectors  computed 
with  spectral  coefficients  at  each  of  the  twelve  levels.  Vectors  were 
computed  with  Eq  (2.3)  on  2\  degree  intervals  in  latitude  and  longitude. 
High  winds  associated  with  the  jet  stream  are  easily  seen  in  Figures  H-7, 

8  and  9  near  tropopause  heights.  The  coefficients  can  be  used  to  compute 
winds  at  any  global  location  where  fallout  or  cloud  transport  predictions 
are  needed. 

Truncation 

Current  NMC  models  use  J  =  30  ,  so  there  are  31  zonal  indices  (A) 
and  32  ordinal  indices  (n)  for  each  Z.  Alternate  truncations  were  examined 
in  this  study  to  determine  how  abbreviation  of  (£,n)  space  affected  wind 
speed  calculations  (4).  Figure  11-14  shows  alternate  rhomboidal  trunca¬ 
tions  of  J  =  27,  24,  21,  18  ana  15.  Wind  vectors  were  computed  at  189 
northern  hemisphere  points  at  the  1000  millibar  level  with  each  truncation, 
noting  the  average  deviation  of  each  component's  magnitude  from  the  cor¬ 
responding  J  =  30  component's  magnitude.  This  is  shown  in  Figures  11-15 
and  11-16.  Deviations  were  largest  when  J  <  24  ,  and  v  component  devia¬ 
tions  were  generally  larger  than  u  component  deviations.  Most  signifi¬ 
cantly,  winds  less  than  1  meter  per  second  (mps)  were  much  more  affected  by 
truncation  than  winds  greater  than  1  mps.  A  truncation  to  J =  24  produces 
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Spectral  Wind  Vector  Fie 


Vector  Field,  400  mb 


Vector  Field,  250  mb 


Vector  Fie 
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Figure  11-10.  Spectral  Wind  Vector  Field,  150  mb 


Figure  11-13.  Spectral  Wind  Vector  Field,  50  mb 
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Figure  11-15.  Truncation  Effect  on  U  Wind  Component 
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Figure  11-16.  Truncation  Effect  on  V  Wind  Component 
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an  average  deviation  of  ±50%  (all  winds) ,  but  only  ±5%  considering  only 
winds  greater  than  1  mps.  Abbreviating  the  summations  would  require 
fewer  calculations  per  component  than  setting  J=30,  but  the  error  intro¬ 
duced,  in  low  velocity  winds  may  cause  large  errors  in  hotline  placement. 
Spectral  wind  calculations  in  this  report  use  J=30. 

Accuracy  of  Spectral  Winds 

The  NMC  routinely  verifies  that  the  spectral  winds  do  represent 
the  radiosonde  data  used  to  derive  the  coefficients  (70) (108) (109) (110) 
(111) .  A  measure  of  accuracy  is  the  root-mean- squared  (rms)  vector  error: 

i  N  1 

Aw  =  ^-7  l  [(u  -  up*  +  (v  -  vr)2l*  (2.4) 

n  L  i=l 

where 

Aw  =  rms  vector  error 

(u,v)  =  spectral  wind  components 

(Uj.,vr)  *  radiosonde  wind  components 

The  rms  vector  error  is  determined  by  NMC  monthly  using  spectral  winds 
and  radiosonde  data.  Table  II-l  shows  the  maximum,  minimum  and  average 
rms  vector  errors,  for  four  levels,  based  on  25  consecutive  months*  data. 
RMS  vector  errors  are  derived  from  data  taken  by  110  North  American 
weather  observation  stations. 

The  vector  errors  do  not  vary  significantly  with  height,  but  wind 
speed  generally  increases  with  height.  Therefore,  the  relative  error, 
ratio  of  vector  error  to  wind  vector  magnitude,  should  decrease  with 
increasing  height.  This  expectation  is  validated  in  Chapter  V,  with  a 
comparison  between  the  spectral  wind  trajectory  and  the  actual  trajectory 
of  a  3  kilometer  high  cloud. 
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TABLE  I I- I 


RMS  Vector  Errors 
(Meters  per  Second) 


Level 

Altitude 

RMS  Vector  Error 

millibars 

kilometers  MSL 

Avg 

Min 

Max 

850 

1.5 

4.8 

4.0 

5.8 

500 

5.6 

4.9 

3.8 

5.9 

250 

10.4 

6.8 

6.0 

7.9 

100 

16.2 

3.7 

2.9 

4.8 

III.  Hotline  Locator  Model 


A  two  step  method  is  used  to  model  variable  wind  effects  with  a 
smearing  fallout  code.  First,  the  hotline  is  located;  then,  activity 
is  smeared  downwind,  along  the  hotline.  This  chapter  describes  how  hot¬ 
line  coordinates  are  computed,  using  spectral  winds.  Essentially,  a 
particle  array  is  tracked  from  the  initial  stabilized  cloud,  through 
spectral  winds,  to  the  ground.  Wind-driven  horizontal  displacements 
of  each  falling  particle  are  computed  and  summed  through  a  discretely 
layered  atmosphere.  The  hotline  is  the  connection  of  particle  landing 
points.  Different  sized  particles  fall  from  different  initial  heights 
with  different  speeds  through  different  winds,  so  hotlines  may  be 
curved. 

Following  is  a  description  of  the  initial  cloud,  atmosphere  dis¬ 
cretization  and  fall-speed  calculation  used  to  locate  the  hotline.  Wind 
shear  is  also  computed  along  particle  trajectories  for  use  in  smearing 
code  estimates  of  cloud  spreading.  Sample  hotlines  are  shown  for  hypo¬ 
thetical  clouds  that  start  above  land-based  missile  sites  in  the  United 
States . 

Initial  Cloud 

Initially,  the  stabilized  cloud  contains  all  radioactive  particles. 
The  cloud  model  is  a  gravity-sorted  distribution  of  particle  size  vs 
height,  obtained  from  correlations  of  DELFIC  cloud  height  data  (11). 
Figure  III-l  shows  how  particle  size  varies  with  height  in  initial 
clouds  produced  by  five  different  nuclear  weapon  yields.  Particle  size 
and  height  in  the  cloud  are  linearly  related;  slopes  and  intercepts 
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are  yield-dependent  polynomials.  Equations  for  particle  height,  slope 
and  intercept  are  in  Appendix  C. 

Assuming  that  particle  starting  heights  represent  altitudes  of 
peaks  in  constant-size  activity  distributions,  the  ground  trace  of 
particles  describes  the  fallout  hotline.  This  assumption  was  validated 
in  (53)  by  comparing  hotlines  from  initial  cloud  correlations  with 
DELFIC  results. 

In  this  report,  initial  clouds  are  the  gravity-sorted  clouds  shown 
in  Figure  III-l.  In  fact,  there  are  many  independent  estimates  of  initial 
cloud  heights,  based  on  numerical  and  empirical  data.  The  distributed 
size-height  correlations  are  used  because  they  replicate  DELFIC 's  initial 
clouds,  and  they  generally  agree  with  observed  data.  Other  initial  clouds 
can  and  have  been  used  for  fallout  calculations  with  analytic  codes. 
Appendix  C  presents  some  initial  cloud  models  and  shows  cloud  height 
sensitivity  to  weapon  yield. 

Atmosphere  Discretization 

The  atmosphere  is  divided  into  layers  of  equal  thickness  beneath 
the  height  of  the  highest  (smallest)  particle  in  the  initial  cloud. 
Particle  trajectory  calculations  showed  that  hotline  location  varied 
with  assumed  number  and  thickness  of  atmosphere  layers  in  the  hotline 
locator  model.  However,  sensitivity  studies  with  3,  6,  12,  24.  36,  and 
48  layer  models  indicate  that  hotline  coordinates  change  very  little 
with  more  than  24  layers. 

Hotline  coordinates  were  computed  with  the  six  different  layer 


models  using  the  same  particle  sizes  in  each  model.  Coordinates  varied 
less  as  atmosphere  discretization  became  finer.  The  deviation  in 


distance  from  the  burst  point  was  used  to  quantify  the  variation  in 
hotline  coordinates  with  number  of  model  layers.  See  the  inset  on 
Figure  III-2.  Average  deviation  is  the  average  of  deviations  of  twelve 
hotline  points.  Figure  III-2  shows  how  the  average  deviation  decreased 
with  increasing  number  of  model  layers.  Hotline  coordinates  were  nearly 
stationary  with  24  or  more  model  layers;  the  average  deviation  due  to 
increasing  the  number  of  model  layers  above  24  was  less  than  4%. 

Atmosphere  properties  are  computed  at  each  level  with  state  equations 
for  a  U.S.  Standard  Atmosphere  (75) (107) .  The  equations  are  summarized  in 
(53) .  Particle  fall  speeds  are  computed  at  each  level  and  averaged  between 
levels  to  obtain  a  fall  speed  and  residence  time  in  each  layer  of  atmo¬ 
sphere.  The  product  of  wind  speed  and  residence  time  gives  displacement 
within  a  layer;  displacements  are  summed  in  all  layers  between  the  par¬ 
ticle’s  initial  height  and  the  earth  surface. 

Wind  speed  at  any  latitude,  longitude  and  height  is  computed  with 
spectral  coefficients  and  Eq  (2.3).  This  hotline  locator  model  linearly 
interpolates  (in  height)  spectral  winds  to  obtain  wind  vector  components 
in  a  layer  of  atmosphere.  The  interpolated  winds  transport,  tailing  par¬ 
ticles  to  new  horizontal  coordinates  at  a  lower  level  in  the  discrete 
atmosphere.  Figure  III-3  shows  the  relations  among  the  discrete  atmo¬ 
sphere,  a  megaton  initial  cloud  and  the  spectral  wind  height". 

Particle  Fall  Rate 

The  terminal  velocity  of  each  particle  must  be  calculated  at 
each  level  to  obtain  the  particle's  residence  time  in  each  discrete 
layer  of  atmosphere.  Terminal  velocities  are  calculated  with  the 
Davles-McDonald  method  (27) (67) (68) .  Davies  empirically  related 
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the  Reynolds  number.  Re,  of  a  moving  sphere  to  CgRe2;  Cq  is  the  particle's 
drag  coefficient.  McDonald  described  an  algorithm  to  compute  the  quantity 
CDRe2  from  particle  size,  particle  density,  atmospheric  density  and  kine¬ 
matic  viscosity.  Using  CQRe2  to  obtain  Re,  terminal  velocity  is  computed 
as  follows: 

V2  =  Re  v/d  (3.1) 

where 

Vz  =  terminal  velocity 

Re  =  Reynolds  number 
v  =  atmospheric  kinematic  viscosity 
d  =  particle  diameter 

The  equations  assume  that  particles  are  spheres  of  constant  density 
falling  through  the  discretely  layered  atmosphere.  Figure  III-4  shows 
the  total  fall  time  of  ten  different  sized  particles  from  heights  up  to 
18  kilometers  above  sea  level  in  a  U.S.  Standard  Atmosphere. 

Wind  Shear 

Vertical  wind  shear  is  a  measure  of  the  wind  field's  dispersive 
effect  on  the  falling  cloud.  Components  of  the  shear  vector  are  defined 
as  follows  (43) (48) : 


=  Au/Az 

(3.2) 

=  Av/Az 

(3.3) 

where 

sx  =  vertical  shear  (westerly) 
s  =  vertical  shear  (southerly) 

Az  =  finite  thickness  of  atmosphere 
(Au,Av)  =  wind  component  variations  over  Az 
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Wind  shear  is  a  characteristic  frequency  with  which  atmospheric  turVu- 

* 

lence  acts  on  the  cloud  (6) (7).  For  nuclear  clouds,  the  turbulence  scale 
is  assumed  to  be  a  fraction  of  the  initial  cloud  thickness.  Fallout 
codes  usually  assume  a  constant  shear  value  for  the  life  of  the  cloud 
(86).  The  constant  shear  is  typically  either  a  guess  or  a  number  cal¬ 
culated  from  the  vertical  profile  that  was  used  to  compute  the  constant, 
'effective*  wind. 

However,  spectral  winds  permit  separate  shear  calculations  for 
each  point  on  the  hotline.  Shear  then  represents  the  net  dispersive 
effect  of  cloud-scale  turbulence  on  individual  particle  size  groups 
in  separate  trajectories.  For  each  hotline  point  (or  particle  size), 
shear  is  computed  as  follows: 


(Sx,Sy) 


(Sx>Sy. 


Ji 


ci 

TFALL 


i 

)2 


(3.4) 


where 


(Sx,Sy)  =  net  (rms)  shear  in  trajectory 
t^  =  residence  time  of  particle  in  layer  i 

TFALL  =  total  fall  time  of  particle  from  initial  cloud  to  ground 


Net  shear  (Sx,Sy)  is  the  rms  value  of  shears  determined  over  discrete 
layers  of  atmosphere,  weighted  by  the  residence  time  of  each  particle 
in  each  layer  of  atmosphere. 


Sample  Hotlines 

Figure  III-5  shows  hotlines  computed  with  a  twenty-four  layer  hot¬ 
line  locator  model  that  tracks  twenty  particles.  The  model  uses  spectral 
coefficients  to  compute  the  winds  that  transport  falling  particles  in  a 
U.S.  Standard  Atmosphere.  Initial  earth  coordinates  of  the  stabilized 
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Curved  Hotlines  with  Variable  Winds  in  Hotline  Locator  Model 


cloud  are  specified  as  colatitude  and  longitude  east  of  Greenwich  in 
radians.  The  hotlines  are  continuous  and  realistically  curved.  Points 
on  the  hotline  are  landing  spots  of  different  particle  sizes;  large 
particles  fall  closer  to  the  burst  point. 

Summary 

A  hotline  locator  model  was  developed  to  compute  hotline  coordinates. 
The  model  uses  spectral  wind  coefficients  to  compute  wind  vectors  and 
shear  along  the  trajectories  of  falling  particles.  Sensitivity  studies 
show  that  computed  hotline  points  converge  to  stationary  locations  as 
the  number  of  atmospheric  layers  increase.  Twenty-four  discrete  layers 
beneath  the  highest  particle  in  the  initial  cloud  give  adequate  conver¬ 
gence.  Wind  shear  is  computed  from  each  particle's  unique  trajectory 
winds.  Sample  hotlines  are  continuous  and  realistically  curved. 

Key  model  assumptions  are  as  follows: 

(1)  U.S.  Standard  Atmosphere 

(2)  Linear  vertical  interpolation  of  winds  between  spectral  heights 

(3)  No  updrafts  or  downdrafts 

(4)  Constant  density  (2.6  grams  per  cubic  centimeter)  spheres 


These  assumptions  were  made  to  generate  the  illustrative  hotlines  shown 
here.  Actually,  any  appropriate  atmosphere,  interpolation  scheme,  verti¬ 
cal  velocity  and  particle  density  can  be  used  to  compute  hotline  location. 


IV.  Variable  Wind  Smearing  Model 

This  chapter  describes  the  second  step  necessary  to  determine  fall¬ 
out  dose  rate  with  variable  winds.  A  fallout  smearing  model  is  presented 
that  determines  ground  coordinates  of  dose  rate  contours  around  curved 
hotlines.  Curved  hotline  coordinates  were  computed  separately  in  the 
first  step  with  the  hotline  locator  model  (Chapter  III),  using  spectral 
winds  for  particle  transport. 

Existing  smearing  codes  use  a  constant  wind  vector  to  transport  a 
falling  cloud  of  radioactive  particles  unidirectional ly  downwind.  The 
one  dimensional  dose  rate  equation  for  smearing  codes  cannot  be  used 
with  variable  winds  because  a  constant  wind  was  assumed  to  obtain  the 
analytic  solution  for  dose  rate  (19).  With  variable  winds,  cloud  acti¬ 
vity  is  smeared  along  a  curved  hotline,  using  a  new  analytic  solution 
to  the  dose  rate  integral. 

Following  sections  describe  the  dose  rate  equation,  its  activity 
arrival  rate  term  and  its  spatial  distribution  term.  A  method  is  devel¬ 
oped  to  compute  dose  rate  at  any  specified  ground  location.  Sample  con¬ 
tours  of  constant  dose  rate  about  curved  hotlines  are  shown. 


Dose  Rate  Equation 

Dose  rate  is  computed  with  an  analytic  solution  to  the  dose  rate 
integral,  Eq  (1.1).  Following  is  the  dose  rato  equation: 


D  =  k  Y  ff  g(ta) 
*^2ir 


(«£  V2  +  a2  VjH  EXP 

A  /  J  A 


(x  v  -  y  vx)2  > 


c2  V2  +  a2  V2  ■» 
x  y  y  x 


in 


(4.1) 


k  =  source  normalization  constant 

Y  =  kilotons  of  weapon  yield 

ff  =  fraction  of  weapon  yield  from  fission 

g(ta)  =  arrival  rate  of  activity  on  ground  (hour-1) 

(X,Y)  =  hotline  coordinates  (miles  from  burst  point) 

(x,y)  =  ground  coordinates  where  dose  rate  is  calculated 
(miles  from  burst  point) 

(V x > Vy)  =  net  wind  components  from  burst  point  to  hotline 
point  (X,Y)  (miles/hour) 

(ax,Oy)  =  deviations  of  cloud  activity  distributions  (miles) 

note:  x  direction  is  positive  eastward 
y  direction  is  positive  northward 


Eq  (4.1)  was  obtained  by  assuming  that  cloud  activity  is  distributed 
normally  in  two  directions,  and  that  the  arrival  rate  function,  g(t), 
can  be  approximated  with  a  two-term  Taylor  expansion  about  arrival  time, 
ta  (19) (53) .  Furthermore,  analytic  integration  of  Eq  (1.1)  required  that 
dose  rate  be  computed  at  points  distant  from  the  burst  point. 

A  constant,  unidirectional  wind  is  simulated  when  Vy=0  and  Vx= 
constant  at  all  hotline  points.  The  dose  rate  equation  then  reduces  to 
the  unidirectional,  constant  wind  form  in  (19);  setting  Vy=0,  the  dose 
rate  equation  becomes: 


6  =  LI  ff  g(ta)  exp  f  -I  1^ 

if  ay  Vx  l  ay  . 


(4.2) 


Activity  Arrival  Rate 

Cloud  activity  is  s/neared  along  the  ground  as  the  falling  cloud 
translates  with  the  winds.  The  activity  arrival  rate  function,  g(t) , 
is  the  functional  rate  of  arrival  of  activity  on  the  ground  (19). 


&(t)  =  -A(r) 


(4.3) 


Assuming  that  a  log-normal  distribution- describes  the  particle  number- 
size  frequency  in  the  cloud,  the  activity-size  distribution  function 
can  be  expressed  as  the  sum  of  two  distributions  (2) (19). 

A(r)  =  A3  +  A2  (4.4) 


A3  =  - — —  EXP 

flU  r 


t  /  m  r  -  <*3  \2 

r*  v  / 


(4.5) 


where 

fv  *  fraction  of  cloud  activity  that  is  volume-distributed  in 
particles 

1-fv  *  fraction  of  cloud  activity  that  is  surface-distributed 
on  particles 
r  =  particle  radius 
8  *  slope  of  log-normal  distribution 

otn  =  aQ  +  n  82  =  logarithm  of  median  radius  in  nth  moment 
distribution 

a0  =  In  r0  =  logarithm  of  median  radius  in  particle  number- 
size  frequency  distribution 

=  time  derivative  of  particle  radius  arriving  on  ground 

A(r)  thus  represents  the  fractionation  or  apportionment  of  cloud  activity 
onto  the  surfaces  (A2)  and  throughout  the  volumes  (A3)  of  the  particles. 
Figure  IV- 1  shows  A3,  ^  and  A(r)  for  typical  values  of  fv  =  0.68  , 

8  *  1.386  and  rQ  *  0.204  (17) (19) (78) .  Figure  IV-2  shows  the  same 

activity  distributions  weighted  by  particle  radius. 

The  time  derivative  of  particle  size  arriving  at  each  hotline  point, 
dr/dt,  is  computed  with  a  Laurent  series  developed  to  compute  the  particle 
size  arriving  from  the  initial  cloud  s,t  specified  times  (24). 


Figure  IV'  1.  Particle  Activity-Size  Distributions  for  Nuclear  Cloud 
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PARTICLE  RADIUS  (MICROMETERS) 


Figure  IV-2.  Particle  Activity-Size  Distributions 
Weighted  with  Particle  Size 


(4.7) 
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where 

=  Laurent  series  constants  (height  dependent) 
t  =  arrival  time  of  particle  size  r  from  a  specified  height 

Each  Laurent  series  constant  (C^  through  Cy)  was  correlated  with  height, 
so  each  constant  can  be  computed  with  a  sixth  degree  polynomial.  The 
polynomials  eliminate  the  need  for  tabular  searches  for  values. 
Appendix  D  contains  Che  polynomial  constants  needed  to  calculate  the 
Laurent  series  constants. 

Since  particles  in  the  initial  cloud  are  gravity-sorted,  different 
particle  Sizes  fall  from  different  initial  heights.  A  different  set  of 
Laurent  series  constants  for  each  particle  size  can  be  computed  with  the 
polynomials.  Thus,  the  distributed  initial  cloud  is  modeled  with  a 
stack  of  pancake  clouds,  one  for  each  particle  size.  The  pancake  clouds 
are  simple,  accurate  approximations  of  vertical  activity  distributions 
that  are  symmetric  about  the  pancake  cloud  height  (12) .  The  pancake 
cloud  computations  of  dr/dt  with  Laurent  series  are  appropriate  only  for 
fallout  particles  with  mass  densities  of  2.6  grams  per  cubic  centimeter. 

For  other  particle  mass  densities,  dr/dt  can  be  computed  with  the 
discrete  particle  arrival  times  generated  by  the  hotline  locator.  Numer¬ 
ical  differentiation  using  r  and  t,  can  produce  any  dr/dt  along  the  hot¬ 
line  (74) .  Figure  IV-3  shows  dr/dt  computed  three  ways  for  a  megaton 
initial  cloud:  (1)  Laurent  series  for  pancake  cloud,  (2)  Laurent 
series  for  a  stack  of  pancake  clouds  at  DELFIC  distributed  heights,  and 
(3)  numerical  derivatives  (central  differences) .  Numerical  values  of 
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dr/dt  were  computed  with  radii  and  arrival  times  from  a  hotline  location 
calculation.  Distributed  clouds  give  essentially  the  same  dr/dt  as  a 
simple  pancake  initial  cloud,  except  at  early  times  (  <  1  hour)  when  large 
particles  are  arriving  from  low  altitudes. 

Spatial  Distribution 

The  spatial  distribution  function  describes  the  lateral  distribution 
of  activity  in  the  cloud. 


Fxv  =  — —  EXP 
/2tt  av 


”  2 


X  v  -  y  V, 


ov 


(4.8) 


_.2  _  —2  \r 2  .  _,2  t/2 

av  ”  ax  Vy  +  °y  Vx 


(4.9) 


Figure  IV-4  illustrates  F^y  over  a  plane,  with  constant  net  wind  components, 
Vx  and  Vy. 

This  two  directional  normal  distribution  defines  the  spread  of  acti¬ 
vity  along  a  line  perpendicular  to  the  net  wind  vector  (8) .  Figure  IV-5 
shows  the  geometry  of  the  net  wind  and  crosswand  lines. 

The  ground  position  (x,y)  is  on  the  net  crosswind  line.  Activity 
is  smeared  simultaneously  onto  the  net  crosswind  lines  as  the  cloud 
moves  downwind  along  the  hotline.  The  spatial  distribution  function, 

FXy,  contains  a  normalized  Gaussian  spread  in  the  net  crosswind  direction. 
Dispersing  clouds  in  the  atmosphere  are  typically  modeled  with  normal 
lateral  distributions.  Data  from  aircraft  sampling  supports  this  Gaus¬ 
sian  model  (29) .  Appendix  F  contains  an  integration  of  FXy  along  a,  net 
crosswind  line,  illustrating  normalization  in  the  crosswind  direction. 
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Figure  IV-4.  Three  Dimensional  Plot  of  Spatial  Distribution  Function 


NORTH 


Net  Wind  and  Hotline  Geometry 


The  deviation,  oy  ,  suggests  an  ellipse-shaped  ground  intercept 
of  the  cloud  with  eccentricity  determined  by  net  wind  components  and 
directional  deviations.  In  fact,  dispersion  theory  and  observations 
of  nuclear  clouds  suggest  an  ellipse-shaped  spatial  distribution.  Air¬ 
craft  sample  data  from  nuclear  clouds  produced  by  Chinese  atmospheric 
nuclear  tests  showed  that  debris  clouds  were  ellipse-shaped  in  horizontal 
directions  at  aircraft  altitudes  (104) . 

Lateral  Deviations 

Standard  deviations  of  the  fallout  cloud  activity  are  specified  in 
the  x  and  y  directions;  x  is  positive  eastward  and  y  is  positive  north¬ 
ward.  The  deviations  are  functions  of  initial  cloud  size,  toroidal 
growth  and  wind  shear  effects.  The  initial  cloud  size  is  determined  by 
weapon  yield.  Toroidal  growth  continues  until  the  cloud's  internal  flow 
stops.  Late  time  cloud  growth  is  driven  by  wind  shear.  Following  are 
equations  for  crosswind  deviation,  Oy  ,  given  by  Pugh  (86)  for  a 
constant  wind,  Vx: 

ay  =  ao  (  1  +  8  )  +  (  t  Sy  ah  )2  (4.10) 


where 

aQ  =  initial  cloud  parameter 


In  0Q  =  0.7  +  -~^  -  3.25  /  [  4.0  +  (  In  Y  +  5.4  )2  ] 


(4.11) 


Y  =  yield  (megatons) 

Ts  =  elapsed  time  for  toroidal  growth  =  t  if  t  <_  3  hours 
Tc  =  effective  e-folding  time  for  activity  arrival  on  the  ground; 
based  on  exponential  g(t)  in  Pugh  (86) 

=  12  (-|j  )  -  2.5  (|§  )2  (hours) 
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(kilofeet) 


HC  =  44.0  +  6.1  In  Y  -  0.205  (In  Y  +  2.42)  | In  Y  +  2.42 { 
ch  =  0.18  HC 
t  =  cloud  arrival  time 

Sy.  =  wind  shear  in  y  direction  (hours" *) 

The  first  term  in  Eq  (4.10)  models  toroidal  growth,  which  is  assumed  to 
stop  at  3  hours  after  the  burst.  The  second  term  represents  growth  due 
to  wind  shear.  Figure  IV-6  shows  how  aQ,  Tc  and  HC  vary  with  weapon  yield. 

With  variable  winds,  downwind  is  not  likely  to  be  in  the  a  direction 
at  all  times.  The  spatial  distribution  function,  FXy,  requires  Oy  and  ox 
to  determine  the  crosswind  spread  of  activity.  Therefore,  ax  is  computed 
with  an  equation  analogous  to  Oy. 

o*  -  crj  (  1  +  8-g  )  +  (  t  Sx  oh  )2  (4.13) 

The  x  and  y  direction  shears  are  computed  from  spectral  winds  using  the 

hotline  locator  model.  Shear  is  determined  along  each  particle's  unique 

trajectory  through  variable  winds. 

In  general  (after  3  hours},  o  and  a  are  proportional  to  t.  This 

X  y 

linear  dependence  is  stronger  than  diffusion  theory  predicts.  Many  models 
of  cloud  dispersion  in  tho  atmosphere  assume  Fickian  diffusion  (20) (96) , 
and  the  analytic  solution  to  the  diffusion  equation  gives  a  Gaussian  form 
with  a  defined  as  follows: 

a  =  (  2  K  t  )  2  (4.14) 

where 

a  =  standard  deviation  of  cloud  distribution 
K  =  coefficient  of  eddy  diffusivity 
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Eddy  diffusivity  coefficients  have  been  estimated  for  nuclear  clouds  (6) 
(64).  With  falling  particles,  the  cloud  dispersion  process  is  affected 
by  vertical  wind  shear,  so  the  diffusivity,  K,  must  be  an  effective  value 
accounting  for  shear,  fall  and  diffusion  (63) (102) (117) . 

Sample  contours  in  this  report  were  computed  with  Pugh's  specifica¬ 
tion  for  Oy,  because  his  formula  was  based  on  nuclear  cloud  data  and 
because  it  contains  explicit  weapon  yield  dependence. 

Dose  Rate  Calculations 

The  variable  wind  smearing  model  determines  fallout  dose  rate  at 
a  point.  Coordinates  of  a  desired  dose  rate  are  found  by  marching  away 
from  hotline  points  (X,Y),  along  net  crosswind  lines,  until  computed 
dose  rate  at  (x,y)  is  the  desired  level.  Figure  IV-7  shows  the  geometry 
for  finding  a  dose  rate  contour;  the  line  connecting  each  point  (x,y) 
is  the  contour.  Figure  IV- 8  shows  0.3  Roentgen/hour  dose  rate  contours 
computed  with  the  curved  hotlines  shown  in  Chapter  III. 

To  compute  dose  rate  at  any  arbitrary  ground  point  (x,y) ,  it  is 
necessary  to  find  the  hotline  point  (X,Y)  that  is  connected  to  (x,y)  on 
a  crosswind  line.  Using  the  fact  that  the  burst  point  (0,0),  hotline 
point  (X,Y)  and  ground  point  (x,y)  form  a  right  angle,  a  method  was 
developed  to  find  (X,Y)  for  any  (x,y) .  Appendix  F  presents  the  technique 
for  finding  the  hotline  point  (X,Y)  associated  with  any  ground  point  (x,y). 
Figure  IV-9  shows  a  three-dimensional  dose  rate  surface  computed  with  this 
technique  at  1000  arbitrary  points  downwind  of  a  megaton  yield  nuclear 
burst . 
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Figure  IV-?,  Dose  Rate  Contour  Geometry 


Summary 

A  method  has  been  developed  to  compute  fallout  dose  rate  contours 
around  curved  hotlines  produced  by  variable  winds.  The  method  uses  an 
analytic  solution  to  the  dose  rate  integral  F,q  (1.1),  and  hotline  coor¬ 
dinates  determined  with  spectral  winds  (Chapter  III).  Activity  is 
smeared  simultaneously  along  net  crosswind  lines  as  the  cloud  moves 


downwind.  Dose  rate  at  any  arbitrary  ground  location  can  be  determined 


V.  Validation 


Ideally,  the  new  fallout  prediction  method  should  be  validated  by 
computing  fallout  dose  rate  contours  for  atmospheric  nuclear  tests,  and 
comparing  calculations  to  actual  data  and  to  other  model  predictions 
(46) (47) (81) .  This  method  requires  NMC  spectral  coefficients  to  compute 
the  wind  vectors  that  transport  falling  particles.  However,  the  NMC  did 
not  generate  and  archive  the  coefficients  before  1980.  Since  the  United 
States  conducted  atmospheric  nuclear  testing  before  1964,  spectral  coef¬ 
ficients  were  not  available  to  validate  the  fallout  transport  method  with 
nuclear  cloud  data. 

The  NMC  Spectral  Model  was  operational  when  the  Mount  St.  Helens 
volcano  erupted  in  1980.  Using  global  wind  data  for  18,  19,  20  May  80, 
(116)  the  NMC  generated  and  provided  spectral  coefficients  to  support 
validation  of  the  new  transport  method  with  ash  cloud  data. 

Ash  cloud  data  for  the  18  May  80  eruption  is  extensive;  radar  obser¬ 
vations,  fallout  collectors,  sampling  airciaft  and  satellites,  provided 
correlative  information  for  fallout  studies.  This  chapter  describes  the 
Mount  St.  Helens  cloud,  the  transport  and  deposition  data,  and  calcula¬ 
tions  of  the  ashfall  hotline  location  and  omars  contours  using  spectral 
winds.  The  excellent  agreement  between  observed  and  calculated  hotline- 
locations  and  ashfall  arrival  times  confirms  that  spectral  winds  can  accu¬ 
rately  determine  trajectories  of  particles  falling  from  a  high  yield 
(me gat on- range)  cloud. 

Further  validation  is  presented  for  a  low  yield  (kiloton  range)  cloud. 
The  Defense  Nuclear  Agency  (DNA)  detonated  600  tons  of  high  explosives  in 
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Now  Mexico  on  26  October  83.  The  test  explosion,  codonamed  DIRECT  COURSE, 
created  a  cloud  that  rose  approximately  3  kilometers  above  ground,  and 
moved  to  the  northwest.  Using  spectral  coefficients,  again  supplied  by 
the  NMC,  the  cloud  trajectory  was  computed  and  compared  to  observations. 
Excellent  agreement  between  calculated  and  observed  trajectories  confirms 
that  spectral  winds  can  accurately  simulate  the  motion  of  low  yield  clouds 
an  the  atmosphere. 

Following  are  the  Mount  St.  Helens  and  DIRECT  COURSE  cloud  analyses. 
Both  cloud  studies  show  that  spectral  wind  transport  is  a  feasible,  accu¬ 
rate  method  for  improving  fallout  modeling  with  variable  winds-  The  Mount 
St.  Helens  analysis  was  extended  to  show  how  the  smearing  model  can  be 
adapted  to  compute  volcanic  ashfall  and  to  back- calculate  the  cloud  par¬ 
ticle  size  distribution  from  ashfall  data. 

Mount  St.  Helens  Ash  Cloud  Analysis 

Cloud  Rise,  Transport  and  Deposition.  The  eruption  of  18  May  80 
produced  the  largest  ash  cloud  of  the  six  major  eruption?  in  1980  (61). 
Initiated  by  a  Richter  5.1  earthquake  at  0832  Pacific  Daylight  Time  (PDT) , 
Mount  St.  Helens  erupted  into  a  Plinian  column,  of  volcanic  gases  and 
particles  that  rose  to  more  than  22  kilometers  above  ground  in  approxi¬ 
mately  10  minutes  (44).  Subsequently,  the  top  part  of  the  ash  column 
expanded  into  a  mushroom- shaped  cloud  that  elongated  and  moved  downwind 
as  ash  particles  fell.  The  volcano  erupted  continuously  for  more  than 
9  hours. 

The  ash  cloud  stabilized  at  different  heights  during  the  nine  hour 
eruption  time.  Table  V-l  shows  the  approximate  altitudes  of  visible 
cloud  top,  bottom  and  center  near  stabilization  time  (90).  Time-varying 
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TABLE  V-l 


Stabilized  Cloud  heights  at  0845  PDT 


Top 

Center 

Bottom 


17.1  kilometers 
12-1 


TABLE  V-2 

Radar-Detected  Cloud  Top  Heights 


Time 

(PDT) 

Height 

0915 

-  0938 

12.0  k 

iloneters 

0938 

-  1100 

12.5 

1100 

-  1158 

14.0 

1158 

-  1255 

13.0 

125S 

-  1330 

12.9 

1330 

-  1500 

13.0 

1500 

-  1605 

13.5 

1605 

15.5 

radar  observations  of  cloud  top  heights  are  shown  in  Table  V-2.  While 
a  visible  cloud  dimension  cannot  be  related  quantitatively  to  a  cloud 
dimension,  the  radar-detected  cloud  top  is  close  to  the  visible  cloud 
center  height.  The  time-weighted  average  of  radar  cloud  top  heights 
is  13.3  kilometers  above  ground  level  (AGL) ,  ranging  from  12.0  to  15.5 
kilometers  during  the  eruption. 

Satellite  data  shows  that  the  ash  cloud  front  travelled  500  kilo¬ 
meters  in  less  than  5  hours.  Figure  V-l  (used  here  with  permission) 
shows  the  time-varying  position  of  the  satellite-detected  cloud’s  lead¬ 
ing  edge;  the  satellite-detected  cloud  was  transported  by  a  high- 
velocity  wind  layer  centered  at  approximately  12  kilometers  above  sea 
level,  A£L.  The  volcano's  mouth  was  approximately  2.9  kilometers  ASL. 
Figure  V-2  shows  the  initial  arrival  times  of  fallout  on  the  ground, 
along  the  downwind  line  of  maximum  ashfall  (the  hotline);  hotline  posi¬ 
tion  was  estimated  from  isomass  contours  in  (90).  Cloud  transport  end 
deposition  are  summarized  in  Figure  V-3  (again,  used  with  permission). 
Ground  ashfall  generally  occurred  to  the  north  of  the  satellite-detected 
cloud,  because  wind  direction  generally  veered  with  increasing  height. 

Figure  V-4  shows  the  elapsed  time  of  eruption,  and  times  for  which 
the  NMC  provided  spectral  coefficients.  Most  of  the  9  hour  eruption  time 
was  within  6  hours  of  OO'JT  19  May  80,  so  cloud  transport  and  fallout  tra¬ 
jectories  were  most-  accurately  replicated  with  persistent  winds  for  that 
time.  Time-varying  winds  were  produced  by  interpolating  spectral  winds 
for  different  times. 

Spectral  Winds.  The  NMC  provided  5  sets  of  spectral  coefficients, 
representing  winds  from  00UT  18  May  89  through  0GUT  20  May  80,  at  12  hour 
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Figure  V-l.  Satellite-Detected  Cloud  Outlines  (90:579) 


DISTANCE  FROM  VOLCANO  (KM) 


Figure  V-2.  Arrival  lime  vs  Downwind  Distance 
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Isomass  Contours  and  Satellite  Cloud  Outlines  (90:588) 


Time  Line  of  Eruption  and  Wind  Coefficients 


intervals.  The  coefficients  were  written  on  magnetic  tape,  formatted  in 
80  column  card  images  (113).  One  set  contains  992  complex  paiTS  for  each 
of  two  wind  vector  components  at  each  of  12  levels.  There  are  23808  com¬ 
plex  pairs  in  each  set  of  coefficients. 

Figures  V-5  through  V-7  show  wind  vectors  computed  with  the  00UT 
19  May  80  coefficients  on  a  2|  degree  latitude/ longitude  grid  over  conti¬ 
nental  United  States.  The  three  levels  were  chosen  to  illustrate  low, 
middle  and  high  altitude  winds  in  the  atmosphere.  High  altitude  tropo¬ 
spheric  winds  are  typically  faster  than  low  altitude  winds  and  the  flow 
is  generally  westerly.  Wind  vector  length  is  the  distance  travelled  by 
a  parcel  of  air  in  3  hours.  The  effect  of  the  jet  stream  is  most  evident 
in  the  250  millibar  wind  field.  Figure  V-6. 

Hotline  and  Isobaric  Traj ectories.  Spectral  coefficients  were  used 
to  compute  tb°  ashfall  hotline  location  and  isobaric  trajectories  at  ten 
heights  above  the  volcano.  An  isobaric  trajectory  is  the  path  travelled 
by  an  air  parcel  at  a  fixed  height  or  pressure  level.  In  1980,  isobaric 
trajectory  predictions  were  used  by  the  United  States  Geological  Survey 
to  forecast  ashfall  hazards  before  the  volcanic  eruption  (73) (99). 

The  hoUJtae  locator  starts  with  an  initial,  stabilized  cloud  and 
computes  the  trajectories  of  spherical  particles  falling  through  vari¬ 
able  spectral  winds  to  the  ground.  The  connection  of  parvicle  landing 
points  on  the  ground  is  the  hotline.  The  initial  cloud  model  was  a 
gravity- sorted  distribution  of  particle  size  vs  height,  described  in 
Chapter  III.  The  size-height  equation  was  derived  from  nuclear  cloud 
calculations,  so  it  is  weapon  yield  dependent.  In  this  application, 
yield  is  just  a  parameter  for  cloud  height  selection.  A  range  of  yields 
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Figure  V-7.  OOUT  19  May  80  Wind  Field  (50  rob) 


was  examined  to  find  a  yield  that  positioned  the  initial  cloud  within 
the  altitude  bounds  detected  visually  and  with  radar,  matching  the 
observed  hotline  location.  Figure  V-8  shows  how  the  computed  hotline 
location  varied  with  different  choices  of  the  yield  parameter  in  the 
OOUT  19  May  80  winds.  The  hotline  from  a  0.5  megaton  initial  cloud 
closely  matched  the  actual  hotline,  although  all  three  of  the  predicted 
hotlines  in  Figure  V-8  have  generally  correct  position  and  curvature. 

The  0.5  megaton  yield  is  not  an  estimate  of  eruption  energy;  yield  is 
a  parameter  used  to  position  the  initial  cloud.  Total  eruption  energy 
has  been  estimated  to  be  much  higher,  approximately  24  megatons  (61:563). 

Furthermore,  computed  arrival  times  of  volcanic  ash  particles  agreed 
with  the  observed  ashfall  onset  times  (38).  Figure  V-9  shows  observed 
and  calculated  arrival  times  on  the  hotline  out  to  1000  kilometers  from 
the  volcano.  Arrival  time  is  an  excellent  indicator  of  the  suitability 
of  spectral  winds  for  fallout  calculations,  because  arrival  time  incor¬ 
porates  all  wind  vectors  computed  along  each  particle's  trajectory. 

Also,  arrival  times  are  computed  during  hotline  location  (Chapter  III), 
and  used  to  define  the  net  wind  vectors  in  the  variable  wind  smearing 
model  (Chapter  IV) .  The  root-mean-square  (rms)  difference  between 
observed  and  computed  net  wind  speeds  is  2.1  meters  per  second.  Figure 
V-9  also  shows  the  arrival  times  computed  with  the  same  initial  cloud, 
using  a  constant  vertical  wind  profile  for  transport  calculations.  The 
constant  vertical  wind  profile  is  a  combination  of  linear  regression 
fits  to  wind  measurements  at  Salem,  Oregon,  at  3  different  times  (21). 

The  constant  wind  profile  generally  overpredicts  arrival  times;  the  rms 
difference  between  observed  and  computed  net  wind  speeds  is  4.3  meters 
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Figure  V-9.  Arrival  Time  Data  and  Calculations 


per  second.  The  spectral  wind  field  produced  net  wind  vectors  that  were 
twice  as  accurate  as  the  net  wind  vectors  computed  with  a  constant  verti¬ 
cal  profile. 

With  spectral  winds  for  different  times,  the  predicted  hotlines 
diverged  significantly.  Figure  V-10  shows  hotlines  predicted  with  spec¬ 
tral  coefficients  for  four  different  times.  A1 3  preceeding  calculations 
were  done  with  one  set  of  spectral  coefficients;  winds  were  spatially, 
not  temporally,  interpolated.  The  same  three  dimensional  wind  field 
persisted  during  cloud  transport.  As  expected,  the  OOUT  19  May  80  winds 
most  accurately  replicated  the  observed  hotline. 

Temporal  variations  were  modeled  by  interpolating  spectral  coef¬ 
ficients  or  wind  vectors  for  adjacent  times  as  particles  fell  from  the 
initial  cloud.  Knowing  when  a  particle  starts  to  fall,  the  temporally 
adjacent  coefficient  sets  can  be  identified  and  interpolated  to  compute 
winds  at  any  time  during  particle  fall.  Figure  V-l]  shows  hotlines 
computed  with  temporally  varying  spectral  winds,  using  three  different 
starting  times:  0845,  1145  and  144S  PDT.  The  three  curves  on  Figure 
V-ll  show  how  temporal  changes  in  trajectory  winds  affected  hotline 
placement.  Within  300  kilometers  of  the  volcano,  temporal  wind  varia¬ 
tions  during  the  eruption  did  not  significantly  affect  ashfall  location. 
However,  at  longer  ranges  and  longer  travel  times,  the  effects  of  tem¬ 
poral  wind  changes  are  apparent.  Ash  that  erupted  in  the  first  six 
hours  curved  southward  sooner  than  ash  that  erupted  later,  causing  the 
ashfall  to  be  smeared  over  wider  areas  at  long  distances.  Temporally 
varying  winds  could  be  used  to  predict  ash  trajectories  if  the  time- 
varying  eruption  rate  data  could  be  estimated  from  radar  measurements 
(44) (45) . 
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Hotlines:  Four  Different  Times 


Hotlines:  Three  Different  Starting  Times  in  Variable 


Small  volcanic  ash  particles,  <1  micrometer  in  diameter  remained 
suspended  in  the  atmosphere  for  a  long  time.  Aircraft  measurements 
showed  that  Mount  St.  Helens  ash  actually  circled  the  earth  in  about 
16  days  (21).  Long  range  transport  was  computed  with  temporally  vary¬ 
ing  spectral  winds  in  constant-altitude  trajectories.  In  this  report, 
constant- altitude  trajectories  are  also  isobaric,  because  the  U.S. 
Standard  Atmosphere  was  used  for  all  atmosphere  state  parameter  calcu¬ 
lations.  Isobaric  trajectories  simulate  the  transport  of  particles 
with  negligible  fall  speeds.  Figures  V-12  and  V-13  show  isobaric  tra¬ 
jectories  near  the  volcano  and  over  a  continental  scale,  computed  with 
temporally  interpolated  spectral  winds,  starting  at  0845  on  18  May  80. 
Two  general  paths  are  evident:  (1)  the  upper  tropospheric  winds  moved 
particles  southeastward  over  Wyoming  and  Kansas  before  turning  to  the 
northeast;  (2)  the  low  level  winds  carried  ash  northeastward  into 
Canada.  The  National  Oceanic  and  Atmospheric  Administration  (NOAA)  pre¬ 
dicted  isobaric  ash  trajectories  using  measured  and  forecasted  wind  data 
(26) (101) .  The  NOAA  predictions  are  shown  in  Figure  V-14.  The  same  two 
general  paths  are  also  evident  in  che  NOAA  predictions. 

Ash  Fallout  Model.  The  dose  rate  equation,  Eq  (4.1),  was  modified 
to  compute  mass  of  volcanic  ash  per  unit  area  on  the  ground.  Three 
changes  were  necessary  to  derive  an  ashfall  equation. 

First,  the  source  normalization  constant,  k,  weapon  yield,  Y,  and 
fission  fraction,  ff,  were  replaced  by  the  mass  lofted  into  the  volcanic 
cloud.  Mass  of  downwind  ash  was  approximately  5  x  1014  grams  (44) (90). 
Estimates  of  total  mass  assumed  rock  mass  densities  of  2.0  to  2.8  grams 
per  cubic  centimeter,  representative  of  the  solid,  unfractured  rock  on 


Isobaric  Trajectories  with  Temporally  Varying  Winds 


Isobaric  Trajectories 


the  mountain  (39).  However,  measurements  of  ash  particles'  mass  density- 
showed  many  particles  with  density  less  than  2.0  grams  per  cubic  centi¬ 
meter  (21).  Cloud  particles  car  have  various  mass  densities,  due  to 
different  rock  types  and  becaui  the  particles  can  contain  air  pockets 
and  fissures  produced  during  solidification,  collision  with  other  parti¬ 
cles  and  breakage  during  eruption  (118).  In  this  report,  2.0  grams  per 
cubic  centimeter  was  used  as  a  representative  value.  In  fact,  hotline 
location  was  not  particularly  sensitive  to  density  variation;  see 
Figure  V-15. 

Second,  the  fractionated  activity-size  frequency  distribution, 

A(r) ,  was  replaced  by  a  particle  mass-size  distribution  function  for 
the  ash  cloud.  A  cloud  distribution  was  measured  by  Carey  and  Sigurdsson 
(21)  using  fallout  samples;  a  cumulative  form  of  their  histogram  is 
shown  in  Figure  V-16,  with  the  DELFIC  default  mass-size  spectrum  for  ref¬ 
erence.  The  measured  distribution  suggests  that  the  Mount  St.  Helens  ash 
cloud  size  distribution  was  multimodal.  A  multimodal  size  distribution 
is  consistent  with  qualitative  and  quantitative  analyses  of  ash  fallout 
(16).  Observed  ash  was  a  mixture  of  two  distinctly  different  rock  types: 
old  rocks  from  previous  eruptions  and  new  magma  (21) (39) (91) .  The  mea¬ 
sured  particle  mass-size  frequency  histogram  was  fit  with  an  analytic 
function  using  the  nine  parameter  function  minimization  technique  des¬ 
cribed  in  Appendix  G.  Following  is  a  trimodal  log-normal  function  that 
fits  the  histogram: 


Figure  V-15.  Hotlines:  Particle  Density  Variations 


CUMULATIVE  PROBABILITY 


Figure  V-16.  Measured  Histogram  and  Nuclear  Particle  Size  Distribution 


SI 


where 


M(r)  =  particle  mass-size  frequency  distribution  function 
fj  =  fraction  of  M(r)  composed  of  the  i-th  distribution 
*  log-slope 

«i  =  ln  (tJ) 

rj  =  median  particle  radius  in  the  i-th  distribution 

See  Table  V-3  for  numerical  values  of  distribution  parameters.  Log-normal 
functions  were  selected  because  explosive  processes  often  produce  log¬ 
normal  particle  size  spectra  (2) (40)  and  because  aircraft  measurements  of 
the  ash  cloud  showed  log-normal  distributions  (36).  Figures  V-17  and  18 
show  the  log-normal  fit  and  the  ash  cloud  histogram.  Three  radius  modes 
occur  near  8,  60  and  250  micrometers. 

The  third  change  to  the  dose  rate  equation  necessary  to  derive  an 
ashfall  equation  is  to  replace  yield-dependent  expressions  for  lateral 
deviations  with  empirical  expressions  for  ash  cloud  spreading.  The 
satellite-detected  ash  cloud  expanded  with  time  as  follows  (21): 

w  =  0.006  t  (5.2) 

where 

w  =  width  of  plume  front  (kilometers) 
t  =  elapsed  time  since  start  of  eruption  (seconds) 


Define 

a  =  (  aj  +  ozy  )S 

(5 . 3) 

and 

cvv 

(a  ,ay)  = - T  0 

(  Vx  +  Sy  ) 

(5.4) 

Assume  that 

w  y  4  a 

(5.5) 

So, 

a  =  0.0015  t 

(5.6) 
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Figure  V-17.  Fit  to  Measured  Size  Distribution 
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CUMULATIVE  PROBABILITY 


The  preceding  three  changes  convert  the  dose  rate  equation  into 


an  ashfall  equation: 


T  M(r)  1 dr/dt 
/2tF  a,. 


EXP 


(5.7) 


where 

A  =  ash  mass/area  on  the  ground  at  (x,y) 

T  =  total  ash  mass  lofted  into  initial  cloud  during  eruption 
M(r)  =  particle  mass-size  frequency  distribution 


Figure  V-19  shows  isomass  contours  computed  with  Eq  (5.7)  using 
hotline  coordinates  of  Figure  V-8.  On  the  hotline,  computed  mass  decreased 
with  increasing  distance  from  the  volcano.  Figure  V-20  shows  computed  ash¬ 
fall  on  the  hotline,  compared  to  actual  measurements.  The  fallout  data 
indicates  a  secondary  maximum  at  approximately  325  kilometers  downwind, 
near  Ritzville,  Washington.  Figure  V-3  also  shows  the  "Ritzville  Bulge". 
The  Ritzville  bulge  was  not  predicted  using  the  aeasured  particle  size 
distribution  (Table  V-3)  in  the  ashfall  equation  Eq  (5.7). 


Particle  Sizes  in  the  Falling  Cloud.  The  measured  grain  size  dis¬ 
tribution  represents  the  size  spectrum  of  fallen,  not  falling,  particles. 
The  falling  cloud  particle  size  spectrum  contained  aggregates  of  smaller 
particles;  many  of  the  larger,  falling  particles  were  broken  during 
impact  and  during  the  mechanical  sieving  processes  used  to  size  the  ash 
samples.  Some  falling  ash  clusters  were  captured  and  sized  without  break¬ 
age,  near  Pullman,  Washington;  data  indicates  that  falling  particles  were 
two  to  five  times  larger  than  fallout  particle  sizes  reported  for  that 
location  (52) (100).  The  unbroken  particle  data  is  consistent  with  terMinal 
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Figure  V-20.  Isoraass  on  Hotline 
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velocity  calculations.  Figure  V-21  shows  where  theoretical  fall  calcu¬ 
lations  predict  that  particles  of  various  sizes  and  densities  would  land. 
Computed  sizes  are  consistently  larger  than  measured  sizes  landing  at 
the  same  place. 

Using  Eq  (5.7),  the  actual  ashfall  data  and  the  function  minimiza¬ 
tion  technique  (Appendix  G) ,  a  particle  mass-size  frequency  spectrum  was 
back-calculated  for  the  falling  cloud.  Specifically,  a  size  distribu¬ 
tion  function,  M'(r),  was  found  that  enables  the  ashfall  equation  to 
closely  match  the  ashfall  data,  including  the  Ritzville  bulge  (90). 

Table  V-4  shows  the  parameters  for  M'(r);  the  function  is  plotted  in 
Figures  V-22  and  V-23.  The  size  spectrum  has  a  mode  near  60  ym  and  par¬ 
ticles  smaller  than  10  ym  constitute  less  than  10%  of  the  falling  cloud 
mass,  compared  to  40%  of  the  measured  fallout  mass.  It  is  interesting 
to  note,  however,  that  the  measured  size  distribution  (Figure  V-17)  also 
contains  a  mode  near  60  ym.  The  presence  of  that  60  ym  secondary  mode 
\n  the  measured  size  distribution  could  be  an  artifact  of  its  parent, 
falling  cloud  distribution,  diminished  by  the  particle  breakage  that 
produced  the  dominant  measured  mode  near  8  ym. 

The  back-calculated  cloud  particle  size  spectrum  was  used  in  Eq  (5.7) 
to  compute  ash  mass  on  and  around  the  hotline.  Figure  V-24  shows  the  hot¬ 
line  mass  compared  to  measured  data;  the  Ritzville  bulge  is  accurately 
replicated.  Isomass  contours  in  Figure  V-25  also  show  the  Ritzville  bulge 

with  appropriate  widening  of  the  1.0  gram  per  square  centimeter  contour 

\ 

around  Ritzville. 

At  distances  beyond  500  kilometers  from  the  volcano,  contours  com¬ 
puted  from  spectral  winds  are  slightly  east  of  the  contours  extrapolated 
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Figure  V-21.  Particle  Size  vs  Distance  from  Volcano 
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CUMULATIVE  PROBABILITY 


Figure  V-23.  Calculated  Falling  Particle  Size  Distribution  (Cumulative) 
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Figure  V-24.  Isomass  on  Hotline  (Falling  Cloud  Size  Spectrum) 
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omass  Contours  (Falling  Cloud  Size  Spectrum) 


from  fallout  and  satellite  data  (see  Figure  V-3).  Ashfall  data  from 
South  Dakota  strongly  suggests  that  the  extrapolated  data  contours  in 
Figure  V-3  are  not  accurately  positioned  (28).  In  fact,  the  computed 
contours  in  Figure  V-25  may  better  represent  the  actual  ashfall  locations 
at  those  long  ranges. 

Summary.  Using  spectral  winds  for  particle  transport,  the  Mount 
St.  Helens  ashfall  hotline  was  accurately  computed,  using  a -gravity- 
sorted  initial  cloud  located  at  approximately  12  kilometers  AGL.  Spectral 
winds  also  produced  long  range  isobaric  trajectories  that  agreed  with  NOAA 
predictions  of  ash  transport  over  the  United  States  and  Canada.  Computed 
fallou..  arrival  times  agreed  with  eyewitness  accounts.  An  ashfall  equa¬ 
tion  was  derived  to  predict  volcanic  ash  mass  at  downwind  points.  The 
measured  fallout  mass-size  frequency  distribution  cannot  produce  the  Ritz- 
ville  bulge  in  fallout  calculations.  A  particle  size  distribution  was 
determined  for  the  falling  cloud,  using  the  ashfall  equation  and  ashfall 
data.  The  computed  falling  cloud  spectrum  simulates  the  Ritzville  bulge 
and  suggests  that  a  secondary  mode  in  the  measured  distribution  is  an 
artifact  of  the  falling  cloud's  mass-size  frequency  distribution.  Figures 
V-26  and  V-27  show  isomass  surfaces  in  three  dimensions,  computed  with 
the  falling  and  fallen  particle  size  distributions,  respectively.  Each 
plot  displays  5000  values  determined  at  surface  positions  that  are  nodes 
of  a  Cartesian  grid,  with  four  kilometer  spacing  laterally.  The  method 
derived  in  Appendix  F  was  used  to  find  the  hotline  points  associated  with 
each  of  the  5000  surface  positions.  The  vertical  axis  is  the  natural  log¬ 
arithm  of  ash  mass,  normalized  to  the  peak  value  near  the  volcano,  and 
truncated  at  1.0%  of  the  peak  value. 
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DIRECT  COURSE  High  Explosive  Cloud  Analysis 


Description  of  Test  Event.  A  six  hundred  ton  sphere  of  high  explo¬ 
sives  was  detonated  on  26  October  83  at  the  White  Sands  Missile  Range. 

The  explosives  were  suspended  from  a  fiberglass  and  steel  tower,  approxi¬ 
mately  fifty  meters  above  ground  level  (AGL)  to  simulate  a  weapon  airburst. 
The  explosion  produced  a  particle  cloud  containing  explosive  by-products, 
dirt  and  tower  debris  that  rose  to  approximately  3000  meters  AGL  as  it  was 
transported  to  the  northwest  by  prevailing  winds.  Climatologically,  winds 
are  westerly,  so  DNA's  debris  (not  radioactive)  fallout  collection  experi¬ 
ments  had  been  positioned  to  the  east  of  ground  zero  before  the  test.  No 
debris  fallout  data  was  obtained.  However,  DNA  tracked  the  cloud  with 
cameras  and  with  a  sampling  aircraft,  so  the  cloud  motion  data  could  be 
used  to  validate  the  spectral  wind  transport  model. 

This  test  is  a  good  validation  exercise  for  spectral  wind  transport. 
The  relatively  low  cloud  simulated  a  cloud  from  a  kiloton  range  nuclear 
burst.  The  coefficients  are  least  accurate  in  the  lower  part  of  the 
atmosphere  due  to  terrain  and  boundary  layer  shear  effects.  Furthermore, 
the  cloud  was  tracked  for  distances  less  than  the  spacing  between  most 
wind  observation  sites  used  to  generate  the  spectral  coefficients.  Thus, 
it  was  theoretically  possible  for  cloud  transport  to  be  affected  by  atmo¬ 
spheric  disturbances  too  small  to  appear  in  the  spectral  model. 

Spectral  Winds.  The  NMC  provided  two  sets  of  spectral  coefficients: 
12UT  26  Oct  83  .nd  00UT  27  Oct  83.  The  explosion  occured  at  18UT 
26  Oct  83.  The  coefficients  were  written  on  magnetic  tape  (114)  in  the 
same  format  as  those  provided  for  the  Mount  St.  Helens  analysis.  Figures 
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V-28  through  V-31  show  spectral  winds  over  the  United  States  computed 
at  the  lowest  four  levels  above  ground  for  GOUT  27  Oct  83.  Easterly 
flow  over  New  Mexico  is  apparently  .caused  by  a  strong  low  pressure  system 
located  over  Northern  Mexico.  Figure  V-32  shows  wind  speed  vs  altitude 
computed  with  coefficients  and  measured  at  the  test  site  (22) .  The 
spectral  winds  lack  the  vertical  resolution  of  radiosonde  data  and  they 
do  not  model  shear  in  the  lowest  layer  of  atmosphere.  However,  because 
the  cloud  transport  was  dominated  by  stronger  winds  at  high  altitudes, 
the  spectral  winds  did  model  the  westward  motion  of  the  cloud. 

Cloud  Motion.  Isobaric  trajectories  are  shown  in  Figures  V-33  and 
V-34  for  both  sets  of  spectral  coefficients.  Figure  V-35  shows  isobaric 
trajectories  using  a  new  coefficient  set,  derived  by  linearly  interpo¬ 
lating  the  12UT  and  OOUT  coefficients  to  represent  winds  at  18UT  26  Oct  83. 
The  aircraft  position  data  illustrates  where  the  sampling  aircraft  was 
located  when  the  pilot  reported  that  he  was  in  any  part  of  the  cloud  or 
stem  (23) .  Table  V-5  shows  the  aircraft  position  data,  including  alti¬ 
tudes.  Spectral  winds  most  accurately  modeled  cloud  transport  near  the 
cloud  top  (500  millibars).  The  spectral  model  clearly  simulated  the 
unusual  easterly  flow  that  dominated  transport  of  the  DIRECT  COURSE  cloud. 

Summary.  The  DIRECT  COURSE  cloud  was  transported  to  the  northwest 
by  climatologically  atypical  winds  (13) (14) .  Spectral  winds  predicted 
isobaric  trajectories  correctly,  to  the  northwest.  The  cloud  trajectory 
and  spectral  wind  trajectory  agree  best  near  the  top  of  the  3  kilometer 
high  cloud.  The  low  altitude  spectral  winds  did  not  model  the  strong 
vertical  shear  in  the  850  to  700  millibar  level. 
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-33.  Isobaric  Trajectories  12UT  26  Oct  83 


Figure  V-34.  Isobaric  Trajectories  OOUT  27  Oct  83 


-35.  Isobaric  Trajectories  18UT  26  Oct  83 


TABLE  V-S 


VI .  Error  Analysis 


Errors  in  the  net  wind  components  Vx  and  Vy  were  shown  in  Chapter  V. 
They  were  determined  by  comparing  the  computed  and  actual  arrival  times 
of  fallout  from  the  Mount  St.  Helens  ash  cloud.  The  computed  rms  error 
of  2.1  meters  per  second  in  the  net  wind  was  even  less  than  the  average 
error  of  5.0  meters  per  second  in  spectral  wind  vectors  (Chapter  II). 

This  chapter  addresses  how  uncertainties  in  the  spectral  winds  contribute 
to  error  in  computed  dose  rate. 

Sources  of  Error 

The  error  in  a  dose  rate  calculation  is  influenced  by  uncertainties 
in  the  independent  variables  of  the  dose  rate  equation,  Eq  (4.1).  Stand¬ 
ard  deviations  of  several  variables  are  listed  in  Table  VI-1;  deviations 
were  based  on  the  ranges  of  published  values.  Table  VI-1  also  shows  best 
estimates  of  the  variables  used  in  this  analysis. 

Error  Propagation  Formula 

An  error  propagation  formula  for  the  dose  rate  equation  was  derived 
by  assuming  that  errors  are  symmetric  about  zero,  and  that  errors  in  any 
two  independent  variables  are  uncorrelated  (10) . 
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TABLE  VI- 1 


Nominal  Values  and  Standard  Deviations 
of  Variables  in  Dose  Rate  Equations 


VARIABLE 

NOMINAL 

0 

REFERENCES 

k, source  normalization 

2350  R  mi2 

0.15  k 

(12) (19) 

constant 

hr  kt 

(41)  (86) 

ff, fission  fraction 

0.5 

0.5  ff 

(19)  (31)  (41) 

a0, logarithm  of  median 

0.203 

2.0 

(12) 

particle  radius 

B3=32> logarithmic 

1.386 

1.0 

(12) 

slopes  of  activity- 

size  distributions 

fv .volume  distributed 

0.68 

0.1  fv 

(12) 

activity  fraction 

(Vx,Vy),net  wind  vec¬ 

computed  with 

5  mps 

(70) (108) (109) 

tor  components 

spectral  winds 

(110) (111) 

(Sx,Sy) ,wind  shear 

computed  with 
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The  error  propagation  formula  Eq  (6.1)  describes  how  errors  may 
propagate  from  independent  variables  into  a  dose  rate  error.  Generally, 
the  variances  are  squares  of  the  deviations  estimated  from  published 
uncertainty  bounds.  Analytic  expressions  for  partial  derivatives  were 
obtained  by  differentiating  the  dose  rate  equation  with  respect  to  each 
independent  variable.  The  dose  rate  equation  depends  on  net  wind  expli¬ 
citly  in  the  spatial  distribution  term  and  implicitly  in  the  arrival  rate 
term. 

Dose  Rate  Error 

A  hypothetical  megaton  burst  was  studied  to  quantify  the  error  in 
dose  rate  caused  by  uncertainties  in  the  independent  variables.  Table 
VI-2  shows  the  relative  and  absolute  magnitudes  of  wind  component  errors 
that  contribute  to  dose  rate  errors.  Values  are  averages  for  ten  ground 
points  located  on  the  200  Roentgen  contour.  The  200  Roentgen  exposure 
contour  was  chosen  because  200  Roentgens  produce  approximately  200  rems 
dose  to  the  surface  of  a  human  body  from  .3  to  3  MeV  gamma  rays.  A  200 
rem  dose  is  the  threshhold  for  radiation  sickness,  based  on  data  from 
Hiroshima  and  Nagasaki  (41).  Average  dose  rate  on  the  contour  is  42.3 
R/hr  at  one  hour  after  burst. 

Spectral  wind  uncertainties  contribute  minimally  to  the  error  in 
computed  dose  rate.  About  4%  of  the  dose  rate  variance  can  be  attri¬ 
buted  to  wind  uncertainty  when  spectral  coefficients  are  used  for  parti¬ 
cle  transport  calculations.  This  analysis  assumed  that  the  deviation  in 
net  wind  components  was  5  meters  per  second.  The  net  wind  vectors 
computed  in  the  Mount  St.  Helens  analysis  diverged  from  actual  data  by 
approximately  2  meters  per  second.  Therefore,  a  spectral  wind  component 
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,  TABLE  VI- 2 


•  Wind  Error  Contributions 


u  =  V 

X 

U  =  Vy 

Dose  Rate  Error 

°u 

(25) 

'•3uJ 

33. 

30. 

Relative  Error 

au 

(25) 

0.80 

0.74 

• 

D 

Fraction  of  Dose 

°u 

,30.2 

'•3uJ 

0.043 

0.039 

Rate  Variance 

_ 

a? 

D 

may  contribute  even  less  than  4%  of  the  dose  rate  variance. 

The  relative  errors  are  consistent  with  similar  error  estimates 
determined  by  W.  Slinn  in  his  assessment  of  meteorological  uncertain¬ 
ties  affecting  the  Nuclear  Regulatory  Commission  Reactor  Safety  Study 
(98).  Slinn  estimates  that  the  relative  error  in  radionuclide  concen¬ 
trations  downwind  of  an  accidental  reactor  release,  due  to  wind  speed 
uncertainties ,  is  approximately  50%  of  the  mean  concentration.  This 
error  analysis  indicated  that  wind-induced  dose  rate  errors  are  80% 
(maximum)  of  the  mean  dose  rate. 
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VII.  Summary,  Conclusions  and  Recommendations 


This  chapter  summarizes  the  analyses,  findings  and  conclusions  of 
this  research.  In  general,  the  study  demonstrated  the  feasibility  of 
using  spectral  coefficients  to  compute  fallout  particle  transport  any¬ 
where  in  the  atmosphere.  Spectral  winds  will  improve  state-of-the-art 
fallout  modeling  by  accounting  for  wind  variability  effects  on  cloud 
trajectories.  This  new  technique  was  validated  by  matching  a  volcanic 
ash  cloud  hotline  and  a  high  explosive  cloud  trajectory.  Wind  errors 
were  actually  quantified  by  comparing  observed  and  calculated  particle 
arrival  times  from  the  volcanic  ash  cloud. 

Particle  Transport  with  Spectral  Winds 

The  NMC  spectral  coefficients  were  used  with  truncated  orthogonal 
polynomials  to  compute  wind  vectors  along  trajectories  of  falling  parti¬ 
cles.  Computed  winds  were  linearly  interpolated  between  the  discrete 
spectral  heights  and  between  sequential  times.  Particle  trajectories 
were  smooth  and  continuous  in  all  dimensions.  Reducing  the  truncation 
limit  to  fewer  than  30  waves  significantly  altered  winds  with  speeds 
less  than  1  meter  per  second.  This  study  used  data  derived  from  radio¬ 
sonde  fits;  forecasted  coefficients  could  be  used  to  predict  fallout 
transport  if  real  time  data  is  not  available  (e.g.  post-attack). 

Variable  Wind  Smearing  Model 

Variable  winds  produced  realistically  curved  hotlines,  so  a  smear¬ 
ing  code  was  developed  that  incorporated  the  two-step  method  to  deter¬ 
mine  dose  rate  at  a  point. 
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Hotline  Location.  A  hotline  locator  model  was  constructed,  using 


spectral  winds  to  transport  an  array  of  trace  particles  falling  from 
the  initial  cloud  to  the  ground.  The  hotline  was  defined  as  the  piece- 
wise  linear  connection  of  particle  landing  points  on  the  ground.  The 
initial  cloud  was  a  correlation  of  DELFIC's  gravity- sorted  particle 
wafer  heights.  Sample  hotlines  were  computed  for  several  locations. 
Wind  shear  was  also  computed  with  spectral  winds  along  each  particle's 
trajectory. 

Dose  Rate  Calculation.  A  fallout  smearing  model  was  developed 
for  use  with  curved  hotline  coordinates  that  were  generated  by  the  hot¬ 
line  locator  and  spectral  winds.  The  model  determines  the  coordinates 
of  any  specified  dose  rate  contour  around  a  hotline.  A  technique  was 
devised  to  compute  dose  rate  at  arbitrary  points  downwind  of  the  burst. 
The  two  dimensional  lateral  spread  of  cloud  activity  is  simulated  with 
a  bi-directional  Gaussian  function,  normalized  in  the  net  crosswind 
direction.  The  dose  rate  equation  for  variable  winds  reduces  to  the 
constant  wind  equation  when  a  constant  unidirectional  wind  vector  is 
assumed  for  cloud  transport. 

Mount  St.  Helens  Analysis 

Ash  cloud  data  from  the  18  May  80  eruption  of  the  Mount  St.  Helens 
volcano  was  used  to  validate  the  new  fallout  transport  method.  With  an 
initial  cloud  that  was  positioned  within  the  observed  heights  of  the 
stabilized  eruption  plume,  the  ashfall  hotline  location  was  computed 
with  spectral  coefficients  provided  by  the  NMC.  The  computed  hotline 
location  accurately  replicated  the  actual  fallout  hotline  and  correctly 
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curved  to  the  southeast  through  Idaho  and  Wyoming.  Hotline  position 
was  not  sensitive  to  particle  density.  Computed  arrival  times  of  vol¬ 
canic  ash  at  points  on  the  ground  were  compared  to  observations,  this 
comparison  provided  an  excellent  quantification  of  the  uncertainty  in 
the  net  wind  vector  computed  with  spectral  winds. 

Hotline  calculations  with  spectral  coefficients  for  four  different 
times  showed  that  particle  transport  was  dominated  by  winds  reported  for 
OOUT  19  May  80.  Temporally  varying  spectral  winds  were  computed  by 
interpolating  spectral  winds.  Those  calculations  showed  that  part  of 
the  plume's  lateral  spread  was  caused  by  time-varying  winds  during  the 
nine  hour  eruption.  Isobaric  trajectory  calculations  with  spectral 
winds  also  agreed  with  published  data,  showing  low  altitude  ash  going 
eastward  over  Canada  and  high  altitude  ash  over  the  United  States. 

A  volcanic  ashfall  equation  was  derived  from  the  dose  rate  equation. 
Mass  per  unit  area  of  volcanic  ash  was  computed  along  the  hotline,  and 
isomass  contours  were  determined  using  a  measured  ash  particle  size 
distribution.  Computed  contours  were  similar  to  fallout  data,  but 
they  did  not  replicate  the  Ritzville  bulge.  Using  the  ashfall  data, 
a  particle  size  distribution  was  back-calculated  that  did  reproduce  the 
mass  bulge  in  isomass  contours  and  on  the  hotline.  The  back- calculated 
distribution  had  a  single  mode  near  60  um;  the  measured  distribution 
was  trimodal,  with  modes  near  8,  60  and  250  ym.  The  60  ym  mode  in  the 
measured  distribution  may  be  an  artifact  of  the  actual  falling  cloud's 
size  spectrum,  before  collision,  impact  and  mechanical  sieving  broke 
the  aggregate,  clustered  particles  into  smaller,  constituent  pieces. 

By  analogy,  particle  breakup  may  also  have  happened  in  nuclear  fallout 
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particle  measurements,  especially  in  Pacific  tests  where  large  amounts 
of  water  were  entrained  into  the  cloud,  facilitating  aggregation  of 
small  particles  into  larger  clusters. 

DIRECT  COURSE  Analysis 

The  actual  cloud  trajectory  was  determined  from  visual  and  air¬ 
craft  tracking  data.  Isobaric  trajectories  were  computed  with  two  sets 
of  spectral  coefficients  supplied  by  the  NMC.  Spectral  winds  correctly 
specified  the  northwestward  cloud  motion.  However,  comparisons  of  actual 
wind  profiles  with  spectral  wind  profiles  demonstrated  that  spectral 
winds,  as  expected,  lack  vertical  resolution  and  shear  between  the  low¬ 
est  spectral  height  and  the  ground.  Accordingly,  the  spectral  winds 
more  accurately  replicated  the  trajectory  of  the  top  half  of  the  high 
explosive  cloud  than  the  bottom  half.  Particle  transport  with  spectral 
winds  should  be  most  accurate  when  the  initial  cloud  is  so  high  that 
the  earth's  boundary  layer  is  a  small  part  of  the  particles'  trajec¬ 
tories  . 

Error  Analysis 

An  error  propagation  formula  was  used  to  quantify  the  effects  of 
wind  (and  other  variables')  uncertainties  on  computed  dose  rate.  Wind 
component  uncertainties  contributed  approximately  4%  of  the  variance  in 
dose  rate,  assuming  ±5  meters  per  second  for  the  net  wind  uncertainty, 
based  on  errors  in  spectral  wind  vectors.  The  Mount  St.  Helens  data 
suggests  that  the  net  wind  uncertainty  is  actually  closer  to  ±2  meters 
per  second,  so  wind  contributions  to  dose  rate  error  may  be  much  less 
than  4%  of  the  dose  rate  variance. 
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Recommendations 


1.  Use  forecast  coefficients  to  determine  how  far  in  advance  tra¬ 
jectories  can  be  accurately  predicted.  This  study  used  coefficients  for 
winds  at  data  observation  times.  The  NMC  can  derive  forecast  coefficients 
using  a  global  circulation  model. 

2.  Modify  DELFIC  to  use  spectral  winds.  However,  realistically 
variable  winds  may  exacerbate  the  numerical  break-up  problem  previously 
observed  with  DELFIC  at  long  ranges.  Wafers  may  become  widely  separated 
in  variable  winds. 

3.  Develop  a  set  of  spectral  coefficients  from  archived  wind  data 
for  six  nuclear  tests.  The  nuclear  tests  should  be  the  same  ones  pre¬ 
viously  used  to  benchmark  DELFIC,  WSEG,  and  DNAF-1  fallout  codes.  (Test 
Event  Codcnames:  Johnie  Boy,  Jangle-S,  Small  Boy,  Koon,  Zuni,  and  Bravo) 
Use  the  coefficients  and  fallout  data  to  further  validate  the  spectral 
wind  transport  method. 

4.  Develop  a  particle  aggregation  model,  accounting  for  the  key 
processes  that  promote  aggregation  in  particle  clouds. 

5.  In  future  volcanic  ash  fallout  experiments,  measure  the  falling 
particle  size  distribution,  using  nondestructive  techniques,  such  as 
cameras  and  laser  spectrometers, to  discriminate  between  falling  and 
fallen  particle  size  spectra. 
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Appendix  A 


Spectral  Coefficients 

Table  A-l  contains  some  spectral  coefficients  for  the  westerly 
wind  component  on  OOUT  16  Jan  82.  There  are  three  complex  pairs  per 
line,  format  (6E13.7).  The  coefficients  are  listed  for  the  first  two 
zonal  wave  numbers  (£  *  0  and  1) .  The  immaginary  parts  are  zero  when 
l  -  0,  because  only  the  real  part  of  the  summation  is  used  to  compute 
wind  speed,  and  the  spherical  harmonics  are  real  when  l  =  0;  see 
Eq  (2.3) . 


Appendix  B 


Derivation  of  Spectral  Wind  Equations 

Wind  vector  components  at  a  fixed  height  in  the  atmosphere  are 
correlated  in  space  with  the  following  expansion: 

J  j  a | +J+1 

(U,V)  =  (u,v)cos<*>  =  l  l  (u£,v£)  P£(sin<{>)  EXP(i£A) 

i,=-J  n=|a| 

Terms  are  defined  in  Chapter  II. 

The  ^-summation  in  Eq  (B.l)  can  be  expressed  as  follows: 


J  -1  0  J 

I  -  l  *  l  *  l 

&=-j  £=-j  &=o  a=i 


Negative  indices  can  be  rewritten. 


-1  |i|+J+l 

l  l  ,  (uJ,V*)  pJCsin<{))  EXP(iU) 

i=-J  n=|a| 


J  n+J+1 

=  l  l  (Un*,V-*0  P-A  Csin4>)  EXP(-UA) 
4=1  n=£ 


J  4+J+l 

*  l  l  etfV)  V 

£=1  n=£ 


Wind  components  are  real  variables,  so  (95:25) (15:689) : 


Cun*'Vn*)  =  C-D  CU*,V*)* 


Compiex  conjugate  is  denoted  by  .  From  (72:495): 


Yn£  =  (-1  )‘*  (?& 


Therefore,  the  negative  il-summation  becomes: 


J  &+J+1 

l  l  a'S.vj)*  (v^)* 

1=1  Tl=l 


J  £+J+l 

I  I  {  C  Rea^vJ)  ]  pj(sin«  (cosU) 
1=1  n=l 


*  [  -i  Im(Uj.Vj)  ]  P^Csinij))  (-i  siiiJU) 


+  [  -i  Im(U*,V*)  ]  F^Csin<f»)  (cosAA) 


+  [  Re(u£,Vj)  ]  P^(simj))  (-i  sinJU)  } 


U  and  V  are  real  numbers;  the  real  part  of  the  summation  with  conju 
gates  is  the  same  without  the  conjugates.  The  real  part  of  Eq  (B.6)  can 


be  expressed  as  follows: 


J  &+J+1 

Re  I  I  (Uj.vj)  Yj 

1=1  n=£ 


Therefore,  the  spectral  wind  equation  becomes  a  summation  over  positive 
indicies.  Spectral  wind  components  are  determined  with  the  following 
expression: 


(U,V)  =  (u,v)  cos<|> 


J+1  J  &+J+1 

=  l  (UJ,VJ)  Pj(sin*)  +  Re  [  2  j  l  (Uj.vj)  P*(sin<|>)  EXP(iU)  ]  (B.9) 

n=0  £=1  n=£ 


Let  A=1  if  £  =  0  ;  A  =  2  otherwise. 


(U,V)  *  (U,V)COS<j> 


J  £+J+l 

Re  [  A  J  (U*,V*)  P*(sin<|>)  EXP(i£X)  (B.10) 
£=0  n=£ 


Associated  Legendre  Polynomials 

p£  are  associated  Legendre  polynomials,  normalized  as  in  (9). 


pJ(x)  = 


C-l)*+n 
2n  n! 


/  (2n+l)  (n-£)l 
y  2  (n+£) ! 


-  HS'+1 
(1-x2)2  d - 


dx*+1 


Cl-x2)1 


(B. 11) 


To  generate  numerical  values,  the  model  uses  the  recursion  relation  (95): 


x  Pj(x)  =eJ+1Pj+1(x)  +  4  Pj_!Cx) 


(B.12) 


where 


(B. 13) 


124 


Belousov  (9)  shows  the  relation  in  the  form: 


2  a!  x  P™  , 
n  n-l 


00  -  bj  p;_2Cx) 


(B. 14) 


where 


and  bm 
n 


/  (2n+l)  (n-m-1)  (n+m-1)  \$ 
\  (2n-3)  (n-m)  (n+m)  / 


Factoring  2  aJJ  from  both  terms  on  the  right  side  of  Eq  (B.14)  gives  the 
following: 


xp;.iW 


> 


(n-m-l) 

(2n-3) 


(n+m-1) 

(2n-l) 


/  n2~m2  \  i 

\  4n2-l  / 


(B.15) 


Let  n  -+  n+1  and  m  -+  SL  . 


X 

I  V— ' 

1 

f  (n-£)  (n+£)  \i  nSL 

[  (2n-l)  (2n+l)  n-lw 

1 

(  (n+1)2  -  m2  \J 
l  4  (n+1)2  -1  / 

(B.16) 


Substituting  Eq  (B.13)  into  Eq  (B.16),  gives  the  following  form  of  Eq  (B.12) 


p£+iw 


x  P*(x)  -  ej  Pj_iCx) 


en+l 


(B.17) 
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Algorithm 


The  algorithm  for  computing  wind  vector  components  uses  double 
precision  calculations  because  the  finite  summation  involves  addition 
and  subtraction  of  992  terms  per  component. 


1. 

Specify 

<J>  and  X  . 

2. 

Compute 

en  • 

3. 

Compute 

P^+^(sin<j>)  from  the  recursion  relation,  Eq  (B.17) 

4. 

Compute 

EXP(i£X)  =  cosilX  +  i  sin£X  . 

5. 

Compute 

(U,V)  with  Eq  (8.101 . 

6. 

Compute 

(u,v)  =  (U,V)/cos<|>  .. 

Appendix  C 


Cloud  Height  Models 


The  height  of  a  stabilized  nuclear  cloud  in  the  atmosphere  depends 
primarily  on  the  weapon  yield,  atmospheric  properties  that  affect  cloud 
buoyancy  and  winds  above  the  burst.  Generally,  cloud  height  increases 
with  increasing  weapon  yield.  Atmospheric  nuclear  weapons  tests  pro¬ 
vided  data  for  empirical  fits  to  cloud  height.  The  most  recent  and 
comprehensive  correlations  of  visible  cloud  top  and  bottom  height  data 
are  simple  power  laws  (80) . 


^top  =  c 

(C.l) 

'bot  ’  a  “b 

(C.2) 

where 

W  is  kilotons  yield 
H  is  meters  above  ground 


a  =  2228  b 
a  =  2661  b 
c  =  3597  d 
c  =  3170  d 
c  =  6474  d 


0.3463 

W  <  4.07 

0.2198 

W  >  4.07 

0.2553 

W  <  2.29 

0.4077 

2.29  <  W  <  19 

0.1650 

W  >  19 

Yield-dependent  heights  are  plotted  in  Figure  C-l. 

By  comparison,  the  WSEG  and  AFIT  codes  use  the  following  equation 
for  height  of  the  radioactive  cloud  center  (19) (86) . 

HC  =  44.0  +  6.1  In  Y  -  0.205  (In  Y  +  2.42)  | In  Y  +  2.42 j  (C.3) 
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□  MSEC  FALLOUT  CLOUO  HEISHTS,  RET.  86 
O  POWER  LAW  FIT  TO  CLOUO  BOTTOM  OATA,  REF.  80 
A  POWER  LAW  FIT  TO  CLOUO  TO?  OATA,  REF  80 
+  CELTIC  AVERASE  WAFER  HEISHTS,  REF.  S3 
-  NUCLEAR  WINTER  CLOUD  TOPS  AND  BOTTOMS,  REF.lOS 


where 

Y  is  megatons  yield 

HC  is  kilofeet  above  ground 

Figure  C-l  shows  that  WSEG’s  equation  puts  the  radioactive  cloud  near 
the  bottom  of  the  visible  cloud.  Early  nuclear  cloud  data  suggested 
that  the  fallout  cloud  was  in  the  lower  part  of  the  visible  cloud  (92) . 
For  yields  above  1  megaton,  WSEG’s  predicted  cloud  height  is  notably 
lower  than  other  models. 

DELFIC  raises  thousands  of  wafers  filled  with  different  sized 
particles  to  various  heights  by  stabilization  time.  Particle  sizes 
are  generally  gravity-sorted  in  the  initial  cloud.  Figure  C-l  shows 
the  average  height  of  wafers  containing  the  DELFIC  default  distribu¬ 
tion's  activity-median  size  (38  ym  radius)  in  stabilized  DELFIC  clouds. 
The  following  equations  were  derived  from  DELFIC  calculations  of  wafer 
heights  in  stabilized  clouds (11) (53) : 

Hp  *  S  Dp  +  B  (C.4) 

where 

Hp  is  average  height  of  DELFIC  wafer  centers  (meters) 

Dp  is  particle  diameter  (ym) 

S  =  EXP  [  1.574  -  0.01197  In  W  +  0.03636  (In  W)2 

-  0.0041  (In  W)  3  +  1.965xl0~4  (In  W)4  ] 

B  =  EXP  [  7,889  +  0.34  In  W  +  0.001226  (In  W)2 

-  0.005227  (In  W) 3  +  4.17xl0-4  (In  W)4  ] 
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Appendix  D 


Polynomials  for  arent  Series  Constants 

Laurent  series  constants  are  used  to  compute  r  and  dr/dt  in  the 
dose  rate  equation. 


r  =  l  CA  t1"6  +  C7  t”i  (D. 1) 

i=l 

The  constants,  C^,  were  derived  by  Colarco  to  enable  smearing  codes  to 
rapidly  determine  particle  size  arriving  on  the  ground  from  a  specified 
height  at  a  specified  time  (24).  There  are  1750  constants,  7  per  alti¬ 
tude,  250  altitudes. 

Each  of  the  seven  constants,  C^,  was  correlated  with  height  (H^) 
in  kilometers  (3).  Seven  polynomial  fits  were  developed:  one  for  each 
constant.  Each  polynomial  for  a  is  6th  degree  in  height,  containing 
seven  constants  (D^j  ,  j  =  1,7). 


Ci  =  Dil  +  Di2  Hk  +  Di3 


Hk 


+  Di4  ^  *  DiS  +  Di6  +  D 


i7  Hk 


(D.2) 


Table  D-l  shows  the  polynomial  constants,  D^j  .  Figures  D-l  and  D-2 
show  r  and  dr/dt  computed  two  ways:  using  tabular  and  using  Eq  (D.2). 
The  figures  show  a  family  of  curves  for  the  same  initial  cloud  heights 
published  by  Bridgman  (19).  Small  magnitude  coefficients  (C^-Cj)  are  not 
reproduced  by  Eq  (D.2)  as  well  as  large  magnitude  coefficients  (C4-C7)  at 
altitudes  £  4  kilometers. 
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TABLE  D-l 


i  i  i  n  r 


?-  0.0  6.0  12.0  18.0  24.0 

%  TIME  (HOURS) 

% 

fcV, 

^  Figure  D-2.  |dr/dt|  vs  Arrival  Time 


L-- 
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Appendix  E 


Crosswind  Integration  of  Spatial  Distribution 

The  two- directional  spatial  distribution  function  defines  how 
activity  is  distributed  in  the  cloud  horizontally. 


Fxy  =  J  f(x,y,t)  dt  CB.l) 

The  integral  was  solved  by  Bridgman  (19)  using  normal  distributions  in 
two  directions  and  a  constant  unidirectional  wind. 


- i -  EXP 

Sir,  Oy  Vx 


V, 


(E.2) 


The  term  (1/VX)  comes  from  the  transformation  of  the  integration  variable 
from  t  to  x,  and  Fy  is  a  normalized  Gaussian  distribution  in  the  crosswind 
(y)  direction. 


(E.3) 


If  the  wind  is  not  constant  and  unidirectional,  a  more  general  result 
is  derived.  With  variable  winds,  the  spatial  distribution  function  is 
expressed  as  follows: 


(E.4) 
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Figure  E-l  illustrates  the  positions  of  wind  vector  components,  the 
ground  point  (x,y)  and  the  crosswind  direction. 

Integrating  Fxy  in  the  crosswind  direction  gives  the  following: 


L  Fxy  ds  =  1  Fxy  ^dx^  dx 


CE.S) 


If  r  is  a  vector  in  the  crosswind  direction. 


r  =  xi . - xt 

V 

y 


CE.6) 


and  ds/dx  is  defined  by  (58) : 


_  =  (  r  •  r 


(E.7) 


The  integral  becomes: 


°0  / 
/  EXP  -1 

-®  /2ir  a„  ^  \ 


(E.8) 


1  ■  •  ft  1 


1  ”  ( ->  ( 


xVy  ‘ 


CE.9J 


On  CC',  y  =  C-Vx/Vy)  x  ,  so  the  integral  in  Eq  (E.9)  becomes: 


-  f  x*  ( vy  +  k)2  ) 

/  EXP  ' - U-  ] 

-co  l  a2  i 


(E.10) 
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y 


The  above  integral  Eq  (E.10)  is  solved  aralytically  (88): 


/2ir  a. 


vx 

Vv  +  - 
>  v 


(E. 11) 


Therefore, 


(E.12) 


This  solution  reduces  to  (1/VX)  or  (1/Vy)  when  a  constant  wind  is  uni¬ 
directional  along  the  x  or  y  axis.  The  constant  wind  result  was  derived 
by  Bridgman  and  Bigelow  (19) . 
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Appendix  F 


Method  to  Compute  Dose  Rate  at  Arbitrary  Coordinates 

To  compute  dose  rate  at  any  point  (x,y)  downwind  of  the  burst,  it 
is  necessary  to  identify  the  hotline  point  (X,Y)  at  the  intersection  of 
the  net  wind  and  crosswind  lines.  The  point  (x,y)  lies  on  the  crosswind 
line.  See  Figure  IV-5.  This  appendix  describes  a  way  to  find  the  hot¬ 
line  point  associated  with  any  ground  point  by  marching  along  the  hotline 
computing  the  angle  between  the  net  wind  and  crosswind  directions,  until 
a  right  angle  is  found.  The  hotline  point  (X,Y)  is  simply  bracketed 
between  two  discrete  hotline  points,  and  then  determined  explicitly,  with 
out  resorting  to  a  numerical  convergence  scheme. 

Figure  F-l  shows  a  hypothetical  hotline,  with  net  wind  and  crosswind 
lines  drawn  for  three  different  hotline  points: 

(X,Y)  is  the  hotline  point  on  a  net  crosswind  line  to  (x,y). 

(Xi,Yi)  is  the  hotline  point  upwind  of  (X,Y). 

(X^+^,Y^+^)  is  the  hotline  point  downwind  of  (X,Y) . 

The  net  wind  vector  from  burst  point  (0,0)  to  (X,Y)  is  OH:  HF  is  the  net 
crosswind  line  that  is  perpendicular  to  OH.  The  point  (x,y)  is  on  HF. 

The  coordinates  (X,Y)  can  be  determined  by  finding  the  point  (X,Y)  that 
makes  HF  normal  to  OH.  Following  are  line  segments: 

OA  =  (  X?  +  Y?  )2  (F.l) 

AF  =  [  (  Xt  -  x  )2  +  (  Yi  -  y  )2  ]2  (F.2) 

OF  =  (  x2  +  y2  )2  (F.3) 


139 


Figure  F-l.  Hotline  Geometry 
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The  line  OF  is  fixed.  The  angle  0^  changes  to  0^+1  when  (X^.Y^)  changes 
t0  CXi+1.Yi+1).  The  angle  9^  is  determined  from  the  cosine  law  for 
triangles : 


COS0.  = 
1 


+  AF2  -  OF2 
•  OA  •  AF 


(F.4) 


To  find  (X,Y) ,  it  is  necessary  to  march  along  the  hotline,  computing 
cos0  at  each  hotline  point.  When  cos0  >0,  0  =  0^+1  and  0^+1  <  ir/2, 
so  (X,Y),  the  vertex  of  the  right  triangle  OH-HF-OF,  is  between  (X^//^) 
and  (xi+i>^+;L).  Assuming  a  piecewise  linear  hotline,  (X,Y)  can  be  deter¬ 
mined  by  finding  AX,  the  distance  from  X^  to  X;  see  Figure  F-2. 


(F.5) 

(F.6) 

(F.7) 


Since  OH  is  perpendicular  to  HF,  the  slopes  of  both  lines  are  related. 


(F.8) 


X2  -  Xx  +  Y2  -  Yy  =  0 


(F.9) 


Substituting  Eq  (F.5)  for  X  and  Y: 

(XA  +  AX)2  -  (XA  +  AX)  x  +  (Ya  +  m  AX)2  -  (Y.  +  m  AX)  y  =  0  (F.10) 
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A  (AX)2  +  B  AX  +  C  =  0 


CF.ll) 


where 

A  =  1  +  m2 

B  =  2X^  +  2mY^-x-my 
C  =  X?  ♦  Y?  -  x  Xi  -  y  Y1 

The  quadratic  Eq  (F.ll)  can  be  solved  explicitly  for  the  two  roots: 

AX+  =  82  -J&  li 

±  2A 

Since  the  hotline  segment  is  represented  with  line  AB,  the  correct  root 
is  the  value  of  AX  that  places  (X-^  +  AX,  Y^  +  AY)  on  AB.  In  Figure  F-3, 
XjL  +  AX+  is  in  segment  AB.  The  other  root  also  produces  a  right  angle  at 
(X,Y) ,  but  it  does  so  by  placing  (X,Y)  outside  the  AB  segment. 

The  above  method  was  verified  by  computing  dose  rates  at  arbitrary 
coordinates  surrounding  a  curved  hotline  from  a  megaton  cloud.  Figure 
IV-9  is  a  three  dimensional  plot  of  downwind  dose  rates  computed  with 
this  technique.  Peak  dose  rates  occur  on  the  hotline  and  values  decrease 
normally  in  crosswind  directions. 


(F. 12) 
(F.13) 
(F. 14) 


(F.15) 
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X=Xt+AX 


Appendix  G 


Function  Minimization  Method  to  Fit  Particle  Size  Distribution  Data 

This  appendix  describes  a  technique  to  fit  particle  size  distribution 
data  points  with  a  trimodal  log-normal  function.  The  technique  requires 
minimization  of  the  sum  of  the  squares  of  differences  between  data  points 
and  functional  values.  The  technique  was  used  twice  in  this  study:  (1)  to 
fit  points  in  the  measured  size  spectrum  published  by  Carey  and  Sigurdsson 
(21) ,  and  (2)  to  fit  size  distribution  data  points  that  were  back-calculated 
with  the  ashfall  equation  and  actual  fallout  data. 

The  measured  size  distribution  was  digitized  from  the  histogram  shown 
in  (21) ,  yielding  twenty  values  of  grain  size  and  corresponding  weight 
fraction.  At  the  top  of  each  histogram  bar,  the  average  radius  was  com¬ 
puted,  and  used  with  the  frequency  for  that  range  of  particle  sizes. 

Twenty  data  points  were  obtained. 

The  data  points  were  fitted  by  minimizing  the  sum  of  twenty  residuals? 
squared  differences  between  data  points  and  values  of  a  trimodal  log-normal 
function.  See  Eq  (5.1).  The  folio  ing  function  was  minimized  within  a 
nine  parameter  hyperspace  constrained  by  upper  and  lower  bounds  on  the  log¬ 
normal  function  parameters. 

20 

F  =  1  [  M(r,).r  -  ]2  (G.l) 

j=l  J  J  J  J 


where 

N(rj)  =  measured  particle  size  distribution  function  for  twenty 
different  particle  radii 
M(rj)  =  trimodal  log-normal  function 
Tj  -  particle  radius 
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The  distribution  functions  M(rp  and  N(rj),  were  weighted  with  rj  to 
facilitate  selection  of  search  limits  on  the  hyperspace.  This  method 
produced  an  excellent  fit  to  the  measured  data;  the  mean  difference 
between  N(rj)*rj  and  M(rj)*rj  was  0.013  for  the  fit  shown  in  Figure  V-17. 
Log-normal  function  parameters  are  in  Table  V-3. 

Isomass  calculations  with  the  ashfall  equation  and  M(rp  generally 
agreed  with  ashfall  data,  but  they  did  not  replicate  the  Ritzville  bulge. 
However,  preliminary  calculations,  using  the  discontinuous  histogram, 
showed  that  multiple  modes  in  the  mass  size  distribution  produced  multi¬ 
ple  maxima  in  downwind  ashfall.  This  observation  suggested  that  a  dif¬ 
ferent  particle  size  spectrum  may  have  produced  the  bulge. 

Using  fallout  data  and  the  ashfall  equation  Eq  (5.7),  particle 
mass-size  frequency  was  back-calculated  at  13  points  along  the  fallout 
centerline.  The  13  points  were  then  fitted  with  the  trimodal  log-normal 
function,  Eq  (5.1),  using  the  same  nine-parameter  constrained  minimiza¬ 
tion  search  technique  used  to  fit  the  measured  distribution  data.  The 
following  function  was  minimized: 


13 

F  =  I  [  M'Crp-rj  -  N'(rj)-rj  ]2  (G.2) 

j=l 

where 

N'(r^)  =  back-calculated  particle  size  frequency  at  discrete  points 
J  along  the  fallout  centerline 
M' (rj)  =  trimodal  log-normal  function 

This  method  produced  an  excellent  fit  to  the  back-calculated  data;  the 
mean  difference  between  N'(rj)*rj  and  M'(rj)*rj  was  0.021.  The  fit  is 
shown  in  Figure  V-22.  Log-normal  function  parameters  are  in  Table  V-4. 
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